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FOREWORD 


This technical report was written in the Applied Mathematics Section, Science 
and Mathematics Research Group, Warfare Analysis Department. Most of the work 
and all of the machine computations were done in 1974 and 1975. However, the 
essential mathematical steps of the present method were develoy.d in 1968 shortly 
after the publication by D.R. Weidman and M.P. Brunner of their model for the 
solution of the same problem. See NWL Technical Note TN-K/21-68. 
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ABSTRACT 


A solution by deterministic methods is described of the problem of computing 
the single-shot kill probability of a point target at a random point from a uniform 
distribution over the interior of an arbitrary ellipse in the plane, given that the distri- 
bution of shots is uncorrelated bivariate normal with respect to a rectangular coordi- 
nate system in the plane, and that the weapon has a cookie-cutter damage function 
with prescribed lethal radius R. This solution has been programmed at NSWC, 
Dahigren Laboratory. The numerical evaluation of a double integral, whose integrand 
contains the so-called elliptic coverage function, is required. Computer results clearly 
show the superiority of this solution over a non-deterministic, Monte Carlo method 
of Weidman and Brunner. 
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1. INTRODUCTION 


The problem treated in this report is as follows. A bomb falls in an uncorre- 
lated bivariate normal distribution, with standard deviations 0, and o, and mean 
point at the origin, with respect to a rectangular Cartesian coordinate system Chk. 
The bomb is assumed to have a known lethal radius R and a cookie-cutter lethality 
function; that is, everything is destroyed in a circle of radius R with the impact 
point as center, and nothing is damaged outside this lethal circle. This lethality or 
damage function has been used in many weapons effectiveness studies. 


An ellipse of semi-axes a and b, a > b, is given, with center at a known point 
(u, v) in the Chk system and the major axis making a known angle ¢ (phi) with the 
heaxis (see Figure 1). A point target is at an unknown random position in the 
interior of the ellipse, governed by a uniform distribution over the area of the 
ellipse. It is required to find the probability P, that the point target is destroyed by 
a single bomb. 


An expression for P, is derived in terms of a double integral, (2-12) or (2-13), 
which must be evaluated numerically. The expression is a function of eight variables, 
u, v, ¢, a, b, R, o,, 9,, of which seven are independent. We could normalize with 
respect to R, or o,, or in other ways, and thus reduce the number of variables by one. 


The variables of integration in (2-12) and (2-13) are r and 6, parametric coordi- 
nates in the ellipse. In each element of area, the function P(h, k, o,, 9,, R) is 
evaluated by NSWC/DL subroutine PKILL, the coordinates (h, k) of the impact 
point of the bomb being expressed in terms of u, v, c (= b/a), ¢, r and @ (2-11). 
This subroutine is based on [2], and computes the probability that a single bomb, 
falling in an elliptical normal distribution of the type described above, wili fall in a 
circle of center (h, k) and radius R. 


In special cases, namely when and only when o, = 0,, we can use NSWC sub- 
routine CIRCV, computing the circular coverage function [5], instead of PKILL. 
CIRCV is at least ten times faster than PKILL. Subroutine CIRCV computes the 
probability P(R, d) that a single bomb, falling in a standard circular normal distri- 
bution, ic. o, =o, = 1, will fall in a circle of radius R with center offset a distance 
d from the mean point of the distribution. Further details on this special case are 
given in Appendix A. 
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This method, referred to hereafter as the DJH method, in contrast to a different 
approach discussed below, and in [1], is entirely deterministic. The accuracy of the 
computed P, is limited only by that of the numerical evaluation of the double 
integral. If an iterative numerical double integration routine is used with successive 
halving of the step sizes in cach direction, P, can be computed to any reasonable 
accuracy level by continuing the iterations until successive iterates agree with each 
other within a preassigned tolerance. A FORTRAN listing of the program is given 
in Appendix F. 


A non-deterministic, Monte Carlo method is used in [1]. A simulated bomb 
impact point (h, k) is determined by the use of random normal numbers. One 
then computes the area common to the ellipse and the lethal circle, i.e., the circle 
of center (h, k) and radius R. The conditional probability of destruction of the point 
target, given that the bomb hits at (h, k) is then given by the ratio of this common 
area to the total area of the ellipse. This follows fror. the uniform distribution of 
the target positions throughout the interior of the ellipse. By replicating this proce- 
dure a large number of times and taking the mean of the conditional probabilities, 
an approximation p is obtained to the overall kill probability P, . 


The method used in [1] for estimating the error in computing P, is to use the 
probable error (P.E.). However, their expression for this quantity is not correct. The 
correct expression is given by 


PE. = 0.6745. / =P — bY (-1) 
n(n — 1) 


where p, (i = 1,<,....n) is a conditional probability as described above, and f is the 
mean of these n conditional probabilities. This mean, 5, is the estimate obtained 
for the kill probability, Py , by the particular simulation involving the n conditional 
probabilities. 


The validity of (1-1) is established in Appendix C. Although this corrected 
formula for the P.E., being smaller than the incorrect formula by a factor I/fn -1 
(see (C6)), puts their method in a better light, the computer results in Section 3 
show, not surprisingly, the superiority of the method (DJH) of this report for 
computing P,. 


See Appendix C for a detailed discussion of the P.E. in the Weidman-Brunner 
(WB) method. 
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2. DESCRIPTION OF DJH METHOD 


The geometry of the problem is shown in Figure 1. A bomb falls in an uncorre- 
lated bivariate normal distribution centered at C, the origin of the rectangular 
Cartesian coordinate system Chk, with standard deviations 0, and o, along the h- and 
k-axes, respectively. In generalo, > 0,. Ifo, = 0,, we have a special case (Appen- 
dix A). The probability density function of the impact point of the bomb is 


h, k) = 
f( eS 


Figure 1. Geometry of the Problem 


It is assumed that the bomb has a cookie-cutter lethality function with known 
lethal radius R (see first paragraph of Section 1). 


A point target is located at a random point from a uniform distribution over 
the area of a given ellipse, centered at a point O with known coordinates (u, v) in 
the Chk system. The semi-axes are a, b, with a > bin general. If a = b, the ellipse 
is a circle, and this is a special case (Appendix A). The major axis of the ellipse is 
inclined at a known angle ¢ to the h-axis, measured positive counterclockwise from 
the h-axis and in the range —- 90 < ¢ < 90 degrees. 


A rectangular coordinate system Oxy is set up with origin at O, the center of 
the ellipse, and x- and y-axes along the major and minor axes, respectively. 


We shall need the equations of transformation from the Oxy to the Chk 
coordinate system. We set up an intermediate system Ox’y’ (Figure 1), with the 
x’- and y‘-axes parallel to the h- and k-axes, respectively. Then the equations of 
translation between the Chk and Ox’y’ systems, and of rotation through the given 


angie ¢ from the Ox’y’ to the Oxy system, are the following: 
h=x’ + u 
k=y t+v } (2-1) 


eee ee ey (2-2, 


y’ = x sin @ + y cos @ 


Combining (2-1) with (2-2), we get the following equations of transformation froin 
the Oxy to the Chk system, in terms of the constants u, v, and ¢. 


h*=u+t+xcos@ — y sing 


k=v+xsing + y cos¢ a) 
The equation of the ellipse (Figure 1) in the Oxy system is 
x? y? 
x + bt 1 (2-4) 


Let b/a = c, with 0 < c < ! except in the special case in which the ellipse is a circle 
(Appendix A), when c = 1. 


We represent points on or inside the ellipse by curvilinear coordinates, r, 9, 
where 


(2-5) 


enn, 
y = cr sin 0 


The curvilinear coordinates, r, @ can be given geometrical significance in the 
following way. Since b = ac, (2-4) is equivalent to 


y? 
| x? += = g? (2-6) 
} ¢? 
| ° and if we put 
; 
| bax 
once (2-7) 
we have 
| B2 +n? =a? (2-8) : 
or, the ellipse corresponds to a circle of radius a in the &n-plane; while for interior 
points in the ellipse, (2-5) and (2-7) give 
§ =r cos | Py 
n=rsin 0 ae . 
or, r and @ are polar coordinates on and inside the auxiliary circle represented by 
(2-8). 


By a standard theorem in analysis (for example, see [3, p. 348] ), the element 
of area dA or dx dy transforms in the following way 


ax dy 
1a(x, or ar 
dA = ie dr dé = dr dé 
air,6 ax ay 
06 4 «63é@ 


cosé@ csin 0 


; dr dé = cr dr dé 
—rsin@ cr cos @ 


(2-10) 


i8(x,y)/(1,6)| being the Jacobian determinant. 
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Subroutine PKILL, or the elliptical coverage routine, (2], does the following 
calculation. Suppose we are given a rectargular coordinate system Chk and a bomb 
falling in an uncorrelated bivariate normal distribution with mean point C and standard 
deviations o, and o, as in the present problem. Suppose we are also given a circle 
with the center at a given point (h, k) and a given radius R. This subroutine computes 
the probability P(h, k, a,» %,, R) that a single bomb falls in the given circle. This 
probability, computed by subroutine PKILL, is also designated as P,,,,, but Py.) 
should not be confused with the overall single shot kill probability P, computed by 
the DJH or the WB method of this report. Suppose that the lethal radius of the bomb, 
with the so-called cookie-cutter damage function, is also R, say 50 feet (this of course 
is the reason for taking this constant R as one of the input parameters in computing 
Ph, k, Ors Ops R)). Then a puint target at (h, k) is destroyed if, and only if, the 
impact point of the bomb is in the given circle centered at (h, k). It is obvious from 
a simple geometrical figure that the point (h, k) is in the lethal circle of the bomb, 
if, and only if, the bomb’s impact point is in the circle of radius R centered at (h, k). 
Thus, the single shot probability of destruction of the point target at (h, k) is given 
by Pch, k, 0,» 9» R) as computed by subroutine PKILL. 

We proceed to set up the double integral whose value, designated as P,, is 
the probability of destruction of the point target at a random point from a uniform 
distribution over the area of the given ellipse. Given any element of area cr dr dé 
of the ellipse (2-6), the probability that the point target lies in this element of area 
is cr dr d@/(nab), the ratio of the area of the element to the entire area of the ellipse, 
since the distribution of the target position is uniform over the ellipse. If the target 
does lie in this element of area, and if a representative point in the element has 
coordinates (h, k) in the Chk system (all points of the element have the same coordi- 
nates, neglecting infinitesimals of higher order), then the conditional probability of 
destruction of the target is P(h, k, o,,0,,R)* cr dr doj(mab). The coordinates of the 
point being (x, y) in the Oxy system, or (r, @) in curvilinear coordinates, the (h, k) 
coordinates of the point, by (2-3) and (2-5) are given by 


ee) (2-11) 


k=v+rcosé@ sin ¢ + cr sin 6 cos @ 


u, Vv, and ¢@ being constants. Hence the conditional probability of destruction of the 
target, given that it lies in the element of area cr dr dé, is 


P(u +r cos @ cos @ — cr sin @ sin ¢, v +r cos @ sin ¢ 
+ cr sin 9 cos ¢, 0,, 0,, R) cr dr d@/(rab), 


4 
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which could be written more briefly as P(h, k, o,° 9,» R) cr dr dé/(wab), with h, k 
given by (2-11). 


The overall kill probability P, is obtained by integrating these conditional 
probabilities over the entire ellipse, since the target must be somewhere in the 
ellipse. The curvilinear coordinate r goes from 0 to a (every point on the ellipse 
having this coordinate r = a), and @ from 6 to 2”. There is no symmetry here, in 
general, because of the arbitrary position of the ellipse in the Chk system (Figure 
1), so that @ must be integrated through all four quadrants, 0 to 27. The result is 


2a 
ao f ff? 2 +8 000 0 con 9 ~ or sino sin 6, v +1 cos 8 sin ¢ 
na? JO 0 


+ cr sin @ cos ¢, a,,0,, R)+ rdr dé (2-12) 


Here we have replaced c/(rab) by c/(7a + ac) or (1/a?). This can be written more 
compactly as 


l 2 fa 
Ro f ff Pots bs 04s 04) RY edd (2-13) 


with h, Kk giver, by (2-11). The P functions on the right are computed by subrou- 
tine PKILL. The double integral as written in (2-12) or (2-13) is now ready for 
evaluation by a numerical double integration subroutine, which will treat the variables 
r and 9 as dummy variables like the variables of integration in any double inte- 
gral, so that their geometrical significance as curvilinear coordinates in the ellipse, or 
polar coordinates in an auxiliary circle, is no longer relevant to the calculation. 


One fact that we can deduce from (2-12) or (2-13), and which is of some 
theoretical significance, and is useful in designing test cases and criteria for deter- 
mining cases in which P, ~ | or P, ~ 0, is that the computed value of the kill 
probability, P,, is simply the space average, over the ellipse of integration, of the 
values of the elliptical coverage function P(h, k, 0, , 0, , R), appearing in the integrand 
in (2-13) and computed by subroutine PKILL, [2]. 
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We can prove this in the following way. Let B denote the average value of 
P(h, k,o,,0,,R). Then 


B= Sf v0. k, 04.04 Rd faa (2-14) 


where each double integration, in the parametric variables r and 6, is over the ellipse 
of integration. Since dA = cr dr dé by (2-10), we then have 


_ j Pe, k, o,, 0,. R) cr dr dé 
p= (2-15) 
f cr dr dé 


Now we move the two factors of c outside their respective integral signs and cancel 
them, and trivially evaluate the integral 


2n a 
f f r dr dé as za?, 
0 40 


so that we have 
p= * fea k R) r dr dé (2-16) 
=— h, k,o,,0,, R)rdr = 
na? Sf 0 Pe 
which is the same as the expression for the kill probability P, in (2-13), so that we 
have 


P=P, (2-17) 


as stated. 


8. NUMERICAL RESULTS 


The basic equation for the computation of the kill probability P, by the DJH 
method is (2-12), or equivalently, (2-13), with h,k given by (2-11). We have to eval- 
uate numerically a double integral, and we have to call the subroutine PKILL, com- 
puting the elliptical coverage function P(h, k, o,, 0,, R), [2], for each point used in 
the integration. 


The double integration routine used is a CDC software routine, using Simpson’s 
rule in each dimension, provided by Control Data Corporation, manufacturer of the 
CDC 6700 machine. This routine allows for the selection of different tolerances for 
the horizontal and vertical integrations. However, the authors of the present report 
did not give systematic study to the selection of these tolerances until several of the 
cases in the testing program required an excessive amount of computing time, more 
than one minute in some cases (page 16). The error analysis which was then worked 
out, as well as a detailed description of the integration method, is given later in this 
section, pages 21-29. Here we give a brief general description of the method, assuming 
that the positive tolerances, denoted as e, and e,, have been appropriately chosen. 


We follow here the notation of (3-1) and Figure 3 (page 22), the integrand 
F(x, y) being unspecified and the integral being regarded as an iterated one. For each 
integration point y ™ s used in the vertical integration (with respect to y), the single 


integral fore, s) dx or f(s) is approximated by Simpson’s rule. 


We give first an overall summary of the two integrations, followed by a few more 
details, but for fuller details - pages 21-29. For a given y = s, where c¢ s< d, the 
exact value of the integral I F(x, 8) dx is denoted by T,, and the Simpson approxima- 
tion by I. The calculation of I, involves a tolerance which is denoted here as e,. A 
function of y, namely I, or I, is thereby defined on [c, d]. This function, in general, 
has discontinuities at values of y for which the number of subintervals of [a, b] 
changes, say from 16 to 32. subintervals (or from 17 to 33 points), but because of the 
smooth and slowly changing integrand here, (2-12), we can assume that the points of 
discontinuity are finite in number and, therefore, the integral I l, dy is approxi- 


mated by Simpson sums [(d — c)/3n] 2 k, 7 as the number of integration points be- 
comes infinite, the multipliers k, being the familiar Simpson factors 1,4,2,4,...,2,4,1]. 


An overall tolerance e, say ¢ * 0.005, is specified. The objective is to use suffi- 
cient points in both the horizontal and the vertical integrations so that we have a high 
degree of confidence, on heuristic grounds at least, that the final approximation | is 
within ¢ of the true value of the double integral; or, [I — Ti < ¢. The tolerance e, 
referred to above, and a tolerance ¢, used in the vertical integration (with respect to i 
y), depend on e in a manner indicated below. é 


For a given y, or 3, sufficient points on [a, b] are taken to give a high degree of 
confidence that the approximation I, is within e, of the exact value I, or Sr, s) dx. 


[1] 
The usual rule in practice is that two successive Simpson sums bos x k, F(x,, 8), 
n t i270 
and bee ¢ k, F(x,, 8) (where n, ™ 2n,, the x,'s, i™ 0, 1,...,.n, are equally spaced 
My jag 2d 2 1 1 2 


x,’s) agree with each other within ¢,. A stricter criterion could be used if desired, for 
example, that of three successive Simpson sums, the largest and smallest must be 
within ¢, of each other. But, once a criterion has been selected, it will be applicable 
for every value of y or s on [c, d], and a function I, is defined on [c, d]} as described 
above. Note that different numbers of points on [a, b], in general, are used un differ- 
ent horizontal lines, i.e., for different values of y or s. This is already implied by the 
remarks above regarding discontinuities in l, at values of y for which the number of 
points on [a, b] changes. 


{ 
y on [a, b], with the x,"8, j=0, 1,....n, also equally spaced at half the mesh length of the 


It is shown that the integral Fe, dy differs by less than (d — c) e, from the 


exact value of the double integral, Ki, dy or I (3-8), But a sufficient number of 
points on [c, d] are taken so that (to a high degree of confidence at least) the approxi- 


mating Simpson sum dos zk, I, (or $f k, I) is within e, of the value of the 


integral fi, dy (3-10), where e,, like ¢,, is one of the tolerances depending on the 


<3 
overall tolerance ¢. Therefore, the final Simpson approximation, ¢ 3 = zk, I. in the ! 
\ 
{ 
\ 


vertical integration, is within (d — c) €, te, of the exact value I of the double 
integral (3-13). 


This sum, (d — c) e, + ¢,, is equated with the overall tolerance ¢ (3-15). The 
tolerances fy and e, may be chosen in any manner such that this equality is satisfied. : 
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The most obvious choice is €, = ¢,, in which case each of these tolerances has the 
value e/({1 + (d — c}] (3-17). 


We now consider in slightly more detail the computation by Simpson’s rule of 
a typical one of the horizontal integrals (i.c., with respect to x) fig F(x, 8) dx, where 


8 is some fixed value of y on [c, d]. If n subintervals of [a, b) are used, or n + 1 
points, n being an even integer, the Simpson’s rule approximation is 


b-—a 


2, k Fx, 9) = —* Fla, 8) + 4F(a + —, s) 


b-—a 


+2Fa+220= 4% 9) de 4Flat(n — 1) 8] + Fb, of 
n ‘ 


The step size (b — a)/n is repeatedly reduced by successive doubling of n until we reach 
a value such that ) = be n! < 9, in case the usual criterion described above is used. 
A similar method is used in the vertical integration, with tolerance e, . 


First we give the numerical values taken by the input parameters o,, o,, R, u, 
v, a, b, and ¢ in the Table E3 in final form as they appear in this technical report, 
pages 59 through 109. Actually, values ure assigned to the ratio b/a rather than to b 
itself. Later in this section, we discuss the testing program involving both the DJH 
and WB methods, and such topics as simple tests for detecting cases in which P, = 0 
or 1 without performing the double integration of (2-12) or (2-13), and the error 
analysis for the method of double integration as used in the DJH program, in the 
general case where P, is not near 0 or 1. 


The standard deviation o, takes the value | in every case because of the normali- 
zation which is used, a takes the values 0.1, 0.5, 1, and 2 on every page; the ratio b/a 
always takes the values 0.2, 0.5, 0.8, and 1; and the angle ¢ takes the values 0°, 45°, 
90°, and —45° on every page. Therefore, the only parameters whose values can vary 
from page to page are R, o,, u, and v (see Table 1), and these are the parameters 
whose values are listed in the appended table. The total number of cases in the main 
table (Table E3) is 59,904, on 51 pages, 59 to 109. 
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Values of u Values of v 
0.2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 ; 
0.2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 
0.2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 
0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 

(See Note) 
0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
] 0,0.5,1,2,2.5,3,3.5,4, 0,0.25,0.50,1 


4.5,5,5.5,6 


Note: For R = 2, 0, = 0.2: v= 0, 0.25, 0.50, 1.50 
But for R = 2, o, = 0.5, 0.8, 1: v = 0, 0.25, 0.50, 1 


see en ri af, 
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The tolerance, or ¢ in the error analysis, is generally 0.005, indicating that, 
although the main table gives values of P, to three decimal digits, the third digit i 
may be incorrect by several units, and we should round off to two decimal digits to 
obtain reliable results. The tolerance was taken as 0.001 in two pages of the table 
(pages 59 and 63) where o, takes its smallest value, 0.2, and R its two smallest 
values, 0.25 and 0.50. But a tolerance of 0,001 entailed the use of too much compu- 
ter time in general. 


At the top of each page are the values of On %% and R for all of the cases on that 
page, but since o, = 1 always, the significant values are those of ¢, and R. The values 


of u, v and ¢ appear at or near the left end of the lines, and those of a and b as column 


headings. : 
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The angle ¢, on every page, takes the four values listed above (~45°, 0°, 45°, 

90°), and similarly a and the ratio b/a take the four values for each which are listed. 

The most common values for u and v are u = O, 0.5, 1, 1.5, 2, and v = 0, 0.25, 0.5, 1. 

The larger values of u (greater than 2.0) are used with the larger values of R, since 

| for these values of R and small to moderate values of u, the values of P, are generally 

near 1. The kill probability P, increases strictly monotonically with R for fixed 

values of the other input parameters, and we always have P, + 1 as R*. On 

| page 107 of the table, where o, = 1, R = 4, the smallest value of P, foru S$ | is 0.941, 

, and for u < 2 it is 0.838, the smallest value on the page. Hence, in order to cover a 

substantial part of the range (0,1) for P, , an additional page (page 109) was computed 
for o, = 1, R = 4, with u taking values up to 6. 


Also, on pages 87-88 of the tables, the values of v are 0, 0.25, 0.5, 1.5, since on 
these pages, where 0, = 0.2, R= 2, the values of P, are generally insensitive to the 
value of v, and we get a wider range of values of P, with these values of v than with 
the usual values 0, 0.25, 0.5, 1. 


Further comments on the selection of values for the input parameters are given 
later in this section, pages 18-20. 


In checking out the DJH method of this report (see (2-12) and (2-13)), 22 cases 
were run on the CDC 6700 computer. These cases are listed in Table 2. Actually, 
one case was inadvertently assigned two case numbers, 14 and 15, so that there were ; 
only 21 different cases, but we follow the original numbering here. Some of these 
cases were checked by three different methods, (a) double integration using subrou- 
tine PKILL (i.e., the general DJH method), (b) double integration using subroutine 
CIRCV instead of PKILL (see Appendix A, page 31), and (c) single integration using 
CIRCV (see (A2) and (A3) in Appendix A). Methods (b) and (c) can be used if, and 
only if o, = 0,, as stated in Appendix A. 


\ 
i 
1 
{ 
} 
4 


The results were uniformly consistent, justifying the conclusion that the DJH 
program computes P, with an absolute error not exceeding a few units in the third 
decimal digit. Although this level of accuracy is adequate for most applications, five 
to six decimal digit accuracy can be attained for P, by a simple change in the PKILL 
subroutine to yield six decimal digit accuracy instead of three as it now stands, and a 
willingness to pay the costs of significantly increased computer time. The time will, 
in general, be more than doubled, in obtaining results correct to six rather than three 
decimal digits. Subroutine PKILL, [2], on the CDC 6700 machine, requires approxi- 
mately 5 milliseconds for three-digit accuracy and 10 milliseconds for six digits. In 
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Table 2. Test Cases for DJH and WB Methods 


14 


addition, there will be a finer tolerance, say « = 107 °, in the error analysis later in this 
section, pages 21-29, in computing to six digits, than the toleran:e of 10-3 in the 
three-digit calculation, and this will entail more calls to the PKILL subroutine. 


In case 16, computed by the three methods mentioned above, with the objective 
of obtaining six-digit agreement, the smallest and largest of the three results for P,. 
were 0.02251899 and 0.02251922 respectively. We note that method (c), the third 
method used, is quite different from the other two in setting up the integration. 


All except cases 8 through 12 were run by the method of [1], which will be 
referred to as the WB method (Weidman and Brunner). These five cases were designed 
to test the DJH method by using a sequence of decreasing ellipses closing down on 
the point (0,4), the smallest one having a major axis of 0.0625. 


The results obtained by the WB method were of variable quality. In some cases 
they were surprisingly good, as in Case 6 in which the DJH and WB methods gave 


’ results of 0.019962 and 0.019922 respectively. An isolated case of good agreement 


could be due to chance, but in several other cases the agreement was quite good. 


But in some other cases the WB method, when applied more than once, generat- 
ing different simulated bomb impact points each time, gave results that differed signif- 
icantly from one another. Thus, in Case 2, which was run three times with the WB 
and once with the deterministic DJH method, the WB method, using 25, 50, and 100 
simulated impact points, gave results of 0.011872, 0.005936 and 0.0217015 reapec- 
tively. Yet the apparently correct P,, as given by the deterministic DJH method, 
was 0.044321. Moreover, the 50% confidence interval for the WB method, computed 
by the method described in Appendix C (see also a brief description in the next para- 
graph of this section), in none of the three cases contained the correct value, 0.044321. 
For the calculation in which 100 impact points were used, the probable error was 
0.008998, and so the 50% confidence interval was 0.0217015 + 0.008998 or 
(0.012704, 0.030700). 


The 50% confidence interval is (p — P.E., p + P.E.), where p is the mean of the n 
conditional probabilities {p,} in a particular simulation, and P.E. is the probable error 
computed by (1-1) or (CS). For a justification of this expression for the 50% confi- 
dence interval, see the latter part of Appendix C, pages 45-46 of this report. As for 
the sample of size n = 100, see the passage from [9] quoted in connection with (C7) in 
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this report. The point involved here is that, although the conditional probabilities P, 


in the WB method cannot be assumed to be normally distributed, the means B of 
samples of size n are approximately normally distributed, by the central limit theorem, 
and the larger the value of n, the better is the approximation to normality. Reference 
(9, p. 111] states that this normality assumption is quite accurate in most cases 
for n> 10. Hence, we can assume that, in this appiication, it is very accurate for 
n = 100, and that the failure of the 50% confidence interval to include the apparently 
correct value P, = 0.044321 cannot be attributed to the fact that the means B in the 
WB method are not from an exactly normal distribution. It should be borne in mind 
that, no matter how large the sample size n may be, the probability that the 50% 
confidence interval will cover the true population mean, P, or p,, is, by definition, 
only 0.50. For higher confidence levels such as 90%, the confidence intervals are 
longer than for the 50% level. 


This last-mentioned calculation (WB method with n = 100) used something of 
the order of 36 seconds on the CDC 6700 computer, since this and = similar run, 
using 100 impact points on a different case, took a total of 73 seconds. If fewer 
simulated impact points are used, the machine time is proportionately shorter. But 
since the P.E. (probable error) is approximately proportional to 1/n, where n is the 
number of simulated bomb impact points, the smaller values of n are associated with 
relatively large P.E.’s. 


The computing time per case by the DJH method, using double integration 
(2-12) or (2-13) and subroutine PKILL, was between 1.0 and 1.5 seconds for those 
cases which were timed, of the 22 cases referred to above, in comparing the DJH and 
WB methods. However, the computing time per case was much larger, in some cases 
of the order of one minute or more, in some of the extreme cases which were com- 
puted in selecting a set of input values for the tables in this technical report. This led 
to some changes in the double integration subroutine used, in an effort to increase 
computing speed while maintaining acceptable accuracy. This point is discussed in 
some detail below. 


Simple tests were included in the program for detecting some, but not all, of the 
cases in which the result will be P, < 0.005 and, therefore, P, will be given as 0.00 
in a two-place table, and for detecting some of the cases resulting in P, > 0.995 (and, 
therefore, Pe = 1.00 correct te two decimal places). We will give a discussion in terms 
of a small but unspecified e, although the actual tolerance used in computing Table 
E3 was e = 0.005, so that this table may be regarded as correct to two decimal 
digits. If the values in a table are subject to errors not exceeding 0.005, the table 
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is conventionally regarded as correct to two decimal digits, although there are excep- 
tions, as for instance when the computed value is 0.824, rounding off to 0.82, with an 
etror of 0.003, making the true value 0.827, which rounds off to 0.83 rather than 0.82. 
More generally, if the errors do not exceed 5 x 10-"°—' or (1/2) x 10-", the table 
is conventionally regarded as correct to n decimal digits, although it would be more 
correct to state that the error, after rounding, does not exceed one unit in the nth 
decimal place. 


For the coordinates of the points (h,,,, k,,,) and (h 
these tests, and further details, see Appendix D. 


one K wax) involved in 


In the first of these tests, for the detection of cases in which P, - e and P, 
will be given as zero to the appropriate number of digits, a point (h, ,, K nin) with 
non-negative coordinates is determined, such that, for every point (h, k) on or inside 
the ellipse, we have ihi 2 Dae 4] > kin and, therefore, P(h, k,o,,0,,R) < Pth,,,, 
Kain? %n? %q? R), since the elliptical c average function decreases with increasing |hi and 
also with increasing |k|. Further, since P, as given by (2-12) or (2-13) is simply an 
average value of P(h, k, o,, o,, R) over the interior of the ellipse, by (2-17), it will be 
true in such a case that P, < P(h_,, Kain’ 9p? 9_» R) with strict inequality, since we 
always have strict inequality at interior points of the ellir.«. Hence, if P(h,,. k.i,5 
0,» 9%» R) < e, requiring only one call to the PKILL subroutine for the selected values 
of h and k, we must have P, <e. This test will fail in some cases in its objective of 
saving computing time, since it is possible to have P, < ¢ < P(h is Kuin> Ty» Fn» RI; 
but it cannot lead to an incorrect decision that P, < e, since, if P(h,.. Kin» Oy» 
o,, R) ¢ ¢, we must have P, < ¢. 

We remark here that, in the preliminary computation of some of the pages of the 
present main table with the coarse tolerance of ¢ = 0.050, the criterion P(h,,.k 4. 

0,» 0,, R) < ¢ was replaced by P( hk... 9,,9,»R)< e/2, but this is “on the safe 
side’’ and certainly not less rigorous than the standazd criterion involving e. The objec- 
tive was to get a smoother decrease in the values of P. with increasing u and v, as 
illustrated by the following example. With « = 0.050, sipgose that P(h, Kins Fn> 
0,» R) = 0.045 and P, = 0.040 if computed by the regular DJH program. If the 
regular ¢ criterion is applied, P, = 0.0 (with accuracy to only one decimal digit if 
e = 0.050). But P(h,,,, k.,,.%,> 9,» R) and the computed P, are both greater than 
e/2 = 0.025, so that P, is not printed as 0.0 under the e/2 criterion. The value 
Px = 0.040, though crobebty not accurate to three digits, may have a substantially 
smaller actual error than the maximum error indicated by the analysis, and the value 
0.040 contributes more than the alternative value 0.0 to a smooth decrease in P, 
with increasing u and v. 
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In the analagous test for P, > 1 — ¢, a point (h,,,, k,,) is found, such that, 
max’ Then 


for every point (h,k) on or inside the ellipse, we have [hi ¢ h,,,, iki ¢ k, 

Py > Phas Kuent %n? [+ R) (by reasoning similar to that in the corresponding step 
in discussing the criterion for P, < ¢). If a call to the PKILL subroutine determines 
that P(h..., k 41. 0,,%» R)> I — c, we must have P, > | —e. Here it is possible 
to have P, > 1 —¢ > P(h,,,.k,.,,9,, 9%» R), in which case the test fails to save 
computing time; but if P(h, .., ko..+0,.%» R)> 1 — e, we must have P, > 1 —e 
and, therefore, the result to be printed in the main table will be P, = 1.00 -~ (to 
the appropriate number of decimal digits) without the necessity of evaluating the 
double integral by the regular DJH program. 


Originally, it was planned to compute a table of 16,384 values of P,, with R 
normalized at | and the remaining input parameters as follows: 


a = 1,2,5,10 o,/0, = 0.1,0.5,0.8,1 
b/a = 0.1,0.2,0.5,1 u = 0,2,5,10 

@ = 45°,0°,45°,90° v = 0,2,4,6 

o, = 0.5,1,2,5 R=! 


If all combinations are taken, the resulting table has 4” or 16,384 values. 


However, it was found from theoretical considerations involving somewhat more 
than the test described above for P, > 0.9995, that such a table would contain no 
values exceeding this number (an unrealistically fine tolerance of e = 0.0005 was under 
consideration at this time). A key fact here is that the ratio of R (#1) to the smallest 
value of 0, (0.5) is only 2. 


Several skeleton tables, each of 128 values, were computed with the object of 
selecting a set of input parameters which would result in a better balanced table in the 
sense of having a greater proportion of values close to !.000, and also a reasonable 
proportion close to 0.000. The following is a typical set resulting in 128 cases, with 
R = I in all cases: 


o, = 0.1,2 u = 0,10 
o,/o, = 01,1 v = 0,4 
a=0.1,8 @ = 45°,90° 


b/a = 0.1,1 R=1 


The short tests described above tor detecting cases in which P, < 0.005 or P, > 0.995 
werc used, with positive results in a substantial proportion of cases in these skeleton 
tables, saving computer time. 


Later it was decided to normalize by taking o, = 1 (rather thar R = 1) with 
o, $9, always. A number of sets of 128 test cases were computed with this normal- 
ization, a typical set being: 


o, =1 u = 2.0,3.0 
0, = 0.2,0.5 v= 1.0,1.5 

a = 0.5,2.0 @ = 0°,45° 
b/a = 0.5,1.0 R = 0.5,1.0 


However, in a number of these test cases, the computing time on the CDC 6700 
machine was unexpectedly long, sometimes more than one minute per case. It was 
concluded that the numerical double integration subroutine was probably using an 
unnecessarily fine mesh in many of these cases, and study was given to methods of 
avoiding this. The PKILL subroutine is a relatively slow one, requiring something of 
the order of 5 milliseconds even for 3-place accuracy on the 6700 machine. The 
CIRCV subroutine (circular coverage function) is an order of magnitude faster, but 
can be used if, and culy if, o, = o,, as has been stated. At 5 milliseconds per case, 
only about 200 elliptical coverage functions per second can be computed to three 
place accuracy. Hence, if the double integration subroutine should use a mesh of 128 
subintervals or 129 points on each dimension of the rectangle of integration, there 
would be 129? or 16,641 calls to the PKILL subroutine, requiring more than 83 
seconds in calls to this subroutine at 200 calls per second. 


The error analysis, resulting in improved speed in computing the double integral 
of (2-12) or (2-13), while maintaining acceptable accuracy, is given later in this section. 


A number of test runs of 128 cases at a time, with o, always having the value |, 
were made. Also, curves were computed and plotted in the hk-plane for P(h, k, ,, 
o,» R) = 0.005, and 0.995 for o, = 1, and various values of o, and R. This facilitates 
the determination of ranges of values of the input variables in the double integral of 
(2-12) which would result in values of P, = 0.000 or 1.000 to three decimal digits, 
by the tests described earlier in this section. On the basis of all these tests, it was 
tentatively decided to compute a table of 16 pages with 1,024 values to a page, mak- 
ing a total of 16,384 (= 2'*) cases, with the follcwing input values: 


ee eer 


o, = 1 (1 value) 
6, = 0.2,0.5,0.8,1.0 (4 values) 
u = 0,0.5,1.0,3.0 (4 values) . 
v = 0,0.25,0.50,1.50 (4 values) 
| = 0°,45°,90°,-45° (4 values) 


a = Q.1,0.5,1.0,2.0 (4 values) 
b/a = 0.2,0.5,0.8,1.0 (4 values) 
R = 0.5,1.0,2.0,4.0 (4 values) 


\ 
! 
| 
Taking all combinations, the total number of cases is 47 = %'4 = 16,384 as stated 
, above. 
} 
if 
I 


Later, it was decided to include the value R = 0.25 on the ground that this 
might be a more realistic value of the lethal radius R than the Iorger values, from the 
practical military point of view. This would make a total of 20 pages of tables, or 
5 x 4° = 20,480 cases. 


When the computing of the pages was started, it was found that some of the 
pages had an excessive number of values near 0 or 1, and additional pages were com- 
puted by modifying and augmenting the set of values of the input variables, especially 
Uu, on some of the pages. The resulting total for Table E3 was 51 pages, containing 
59,904 values of P, , as stated earlier in this section, page 11. 


A tolerance « = 0.005 was used on all pages except 59 and 63, on which the 
tolerance was 0.001 (see page 12 for further comments on this). The computing 
time per case varied widely. being 0.10 second or less on some pages, and 0.85 
second or more for other pages. The overall average, from available recorded times, 
was about 0.4 second per case. 


In using the DJH program, based on (2-12), arbitrary positive values can be 
given to o,, 0,, a, b, and R, and arbitrary real values to u, v, and |¢| < 90°. There 
is nothing in the analysis which rules out the case o, > o,. If uor v or both are 
negative, the corresponding positive yalue or values should be taken, with appropriate 
change in 9; for example, u = --4, v = 2, @ = 30° is equivalent to u = 4, v = 2, ¢ = —30°, 
as is clear if a figure is drawn, the elliptical coverage function being symmetric with 
respect to both axes. 
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But if the table in this report is used, with or without interpolation, it must 
be noted that, in the table, 0, = | and, > o,. Hence, it may be necessary to nor- 
malize the input variables (except 9) by making o, = 1, and if 0, >0, originally, the 
roles of h and k must be interchanged, If we have the case o,* 5,0, ™2,u=3,v = 2, 
a= 2, b = 1.5, 6 = 30°, R = 3, we simply divide all distances by 5 and obtain as an 
equivalent case o, = 1. 0, = 0.4, u = 0.6, v = 0.4, a = 0.4, b = 0.3, ¢ = 30° (no change 
in 9), R = 0.6. 


Now, suppose we have a case with o, > o,, say o, = 2,0, =5,u=3,v = 2, 
a= 2, b= 1.5,¢ = 30°, R= 3. We first divide all distances by the larger o, o, in this 
case, and get the equivalent case o, = 0.4, 0, = 1,u = 0.6, v = 0.4, a = 0.4, b = 0.3, 
¢@ = 30°, R = 0.6. Then, reversing the roles of h and k, we will also have to reverse 
those of u and v and make an appropriate change in ¢, obtaining the equivalent case 
o, = 1,0, 70.4, u = 0.4, v = 0.6, a = 0.4, b = 0.3, = 60°, R= 0.6. This should be 
clear from a study of Figure 2, bearing in mind the symmetry of the elliptical coverage 
function with respect to both axes. 
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Figure 2. Interchange of o,, 9, in Example 


Here we discuss the numerical evaluation of the double integral (Figure 3) 


~ af rd 
T= f i F(x. 1 ax ay (3-1) 
c a 


by a method which has been chosen as an efficient one, using two tolerances Ey and 
¢, and an overall tolerance € 
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Figure 3. Area of Integration for Double Integra! 


| 
: 
i} 
| 
} 
These three e’s may be given appropriate positive values depending on the accuracy 
\ required; for example, the overall tolerance e, on which €, and e, depend in the pre- 
f sent analysis (see (3-15)), is sometimes given a value of 0.01 in crude calculations, and 
sometimes a value of 0.00001 if much more accurate results are desired. 

In (2-12) and (2-13), the variables of integration are r and 0. We leave open the 
possibility of identifying r with x and 6 with y, or the reverse. 


We regard the integral as an iterated one, 


-_ rat pb d 
i= f | f F(x, y) wore f fly) dy (3-2) 
c a c i 
| 
} 


b 
For a fixed y, the integral f(y) = J F(x, y) dx is approximated by Simpson’s rule, 


but the number of subintervals for two different ordinates or values of y is not 
necessarily the same. Hence, we do not have a two-dimensional grid of the less 
efficient type sometimes considered in the literature, for example in [7, p. 186, ( 
problem 54}. 


For y = s for arbitrary s on (c, d), let 


I, = exact value of f(y) = f(s) - fF. s) dx (3-3) _ 
A 
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I, = Simpson's rule approximation for T, (3-4) 


In general, a bar will indicate an exact value, as in. 1, in this case indicates an 
approximation. 


I=, +8, (3-5) 


where E, is the error committed by using the Simpson's rule approximation with the 
chosen number of subintervals of [a, b). 


For a given s, i.e., a given value of y, let a sufficiently large number of subinter- 
vals of [a, b] be chosen so that 


1, -Ui iEBi<e, (3-6) 
where ¢, is one of our fundamental tolerances. 


We are really comparing successive approximations, say (I,), and (1,),,, with 
each other, where n or 2n denotes the number of subintervals into which [a, b] is 
divided, rather than comparing any approximation with the exact value Ll; since in 
practice we do not know the value of the latter. We do not in this analysis derive 
rigorous error bounds for integration by Simpson’s rule. Numerical analysts have 
developed several such rigorous bounds, using the fourth derivative of the integrand 
[7, p. 75, (3.5.23)], or evaluating a line integral in the complex plane, or using other 
methods. The calculation of any of these rigorous bounds would require « relatively 
large amount of machine time, and, moreover, such bounds are often ultra-conserv- 
ative, i.e., much larger than the actual error committed. It is felt by the present 
authors that sufficient accuracy for practical computations can be attained by develop- 
ing simple tests involving convergence of successive Simpson approximations to a 
common limit. Such tests can give a high degree of confidence in the reliability of 
the results, though not absolute certainty. We may decide that (3-6) is satisfied if 
two successive iterates agree with each other within ¢,; or, if we are more conserv- 
ative, within ¢,/2; or, of three successive iterates, say CI)» 12, and Ud... the 
largest and smallest are within e, of each other. We assume here that we have some 
such criterion, but we do not specify at this stage what it must be. We state in (3-17) 
below, how the tolerances are handled in the present NSWC program, but the program 
is flexible in that changes in the tolerances are easily made. 
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Now by definition, 


[Mor ftont 


= exact value of double integral (3-7) 


Also, I, is a well-defined function of s or y once a numerical value of ¢, has been 
assigned and a criterion for the satisfaction of (3-6) has been chosen as discussed 
above. 


We assume also that fi, dy or £, dy exists. I, or SOF, y) dx (3-3), the 


exact value of the integral from a to b, is a continuous function of y, by the conti- 
nuity of the integrand of (2-12) as a function of two variables. A consideration of 
the manner of forming the Simpson sums shows that the approximation I, is also a 
continuous function of y in any interval on the y-axis in which the number of sub- 


b 
divisions of {a, b] in forming the Simpson approximation to f F(x, y) dx does not 
change. But l, does have a discontinuity, in general, at a pent where we switch 
from 16 hig 32 subintervals, for example, in approximating I F(x, y) dx. Suppose 
a I F(x, y,) dx uses 16 subintervals for some y = y,. but I F(x, y) dx uses 


32 subintervals for some neighboring values of y, say for all y in the open interval 
(Yy>¥q + k) for some positive k. Then oe as a function of y, will in general have a 
discontinuity for y = y,. Nevertheless, the integrand in (2-12) is such a smooth 
slowly changing function of its arguments that we can be confident on heuristic 


grounds that the integrand in i I, dy has a small finite number of discontinuities. 
But the Riemann integral fA, dy exists if, and only if, the set of discontinuities is 
of Lesbesgue measure zero [8, p. 230). This includes the present case, in on the 
discontinuities are assumed to be finite in number. Hence, the integral I I, dy or 
d 
Je I, dy exists. It is shown below in (3-8) that it differs by less than (d — c)e, 
from the exact value of the double integral, k , dy or I, i.e., it approximates T 
arbitrarily closely for small «,. Also, the Simpson sums vs z k, is which are 
used in approximating the integral with respect to y, the k,’s being the Simpson 
multipliers | ,4,2,....4,1, converge to the value of the integral hk I dy as the number 
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of subdivisions of [c, d] becomes infinite (see (3-10)). It is not difficult to prove that 
the Simpson sums, assuming that they converge to the true value of the integral in the 
integration of a continuous function over a finite interval, also converge to the value 
of the integral if the function has a finite number of discontinuities with a finite jump 
at each such point. 


Assuming then that f I, dy exists, and assuming that (3-6) holds everywhere on 
(c, d], ie., that we take a sufficient number of subintervals of [a, b] so that (3-6) 
holds for all y on [c, d], we have 


fw fi 
sf -4| < fi av=a-o, (3-8) 


fr - fi, ay| <(d-c)e, (3-9) 
fi dy = T 


Now, suppose that we take a sufficiently large number of subintervals on [c, d] 
8o that 


where the sum uses the I,’s for the finite number of points taken on (c,d) (=n+1 
points where n is the number of subintervals of [c, d]). Here, € is another funda- 
reine tolerance, like e,. As in the case of e,, we can require two successive iterates 


,-Wey| 


or 


<e, (3-10) 


is Dk I, ,% agree with each other within e, , or within e ,/2, or choose some other 
criterion. We do not know the value of the integral fi, I_ dy, so we are forced to adopt 
some arbitrary criterion for the satisfaction of (3-10). 
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With (3-6) and (3-10) satisfied, we now have 


d-c¢ d-c : 

[rSStennl = fiers fr-Seeny] ea | 

d-—c d-—c { 

[r- — Ek, | < |t- fi, as| + | fice S22, | (3-12) | 
T Laake <(d—c)e, + 3-13 
ar A, (d—c)e, €, (3-13) 


Now suppose that ¢ is some overall tolerance, and we wish to choose ¢, and ¢, 
in such a way that 


d— 
[r - <r k, I | <e; (3-14) 
3n 
that is, we wish the computed Simpson sum approximating the double integral to be 
correct within e. From (3-13) and (3-14), we equate e with (d — c)e, +¢,, 80 that 
we have 


r-"—* =, | <(d-—c)e, te, me (3-15) 


Equation (3-15) clearly allows some latitude in deciding on the individual values 
of e, and e,. The present NSWC/DL program takes ¢, = ¢,, from which follows 


{1 + (d -—c)] eg ze (3-16) 


| 


-_- . = é€ 
fo 1 T+ 6) 


(3-17) 


If we wish to experiment with various values of the ratio k = ¢, /e,, a simple generaliza- 
tion would be, from (3-15), 


(d-—c)e, +ke, me (3-18) | 
‘ | 

i 3-19; 
“ok +0) ay 
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+ SUN 


and 


ke 


eke =—. 
° k+(d-0c) 


(3-19a) 


where k is any positive real number, less than, equal to, or greater than 1. The case 
; k = 1 reduces to (3-17). It is possible that computational experience would show 
{ that some value such as k = 2 or 1/2 is generally more efficient than k = 1. The 
determination of such an optimal value, or approximately optimal value, of k, would 
\ be purely an empirical matter. We could try various values of k fora fixed value of e, 
and regard as optimal that value of k which resulted in the smallest expenditure of 
machine time while maintaining the specified accuracy. 


A simple modification of this error analysis is necessary in case the constant 
factor 1/(na?) outside the integral signs in (2-12) and (2-13) is not included as part 
of the integrand at every step of the analysis, but only once at the end in computing 
| ; the true value of P,. This is the case in the DJH program. 


Let the integrand in (2-12) and (2-13), not including the factor 1/(ra?), be | 


t 

\ 

denoted by F*(@, r), and, in general, let the asterisk or “star” denote this omission 
of the constant factor. Then (2-13) can obviously be w-itten in either of the forms 


al pas 
P, = I f F(@, 1) do| dr ! 
i a 2n 

= — F%@, 1) d0| 3-20 
oh lh're.9 | a oy | 


where 


F*(6, r) = P(h, k,o,,0,, R)r (3-21) 


n 


F(6, r) = oe F*(@, r) (3-22) 
a2 


and h, k are given by (2-11). | 
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I*= f F* de dr (3-23) 
and 
PL, sI= i I* (3-24) j 
. na? 


Now, following the foregoing analysis but using the “‘star” in the sense which 
has been described, we evaluate the double integral I* to an accuracy such that 


d-c 
3n 


fe - 


zk, Ip] <ee (3-25) 


(with e* to be determined in terms of the overall tolerance ¢), with 


jess oo eK 3-26 
fp ™e l+@-o (3-26) 


as in (3-17). We will then have 


ot iS eet te 


d-—c 
| T- zk, 1,| = 


! 3n na? na? 3n me 
1 d-c 1 
= ao eeu s au g* = 
, 7 T a eh | <me (3-27) 
But we wish to have 
d-c 
: r erga =k, | <e, (3-28) 
and so we put e*/(xa?) = e or 
e* = na? + (3-29) 
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and therefore, from (3-26), 


DI ae So 


td 
€ e 
es = e@ = ————___ = a? 


1 T+a_o) “14+@—0) 


The tolerances are handled in this way in the DJH program. The constant e is set at 
0.005 in the present program, but can easily be changed to some other desired value 
such as 0.001 or 0.05. 
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SPECIAL CASES 


Two special cases are worthy of mention here: (a) 0, =9,;(b) a=b, ile, 
the ellipse is a circle. 


(a) We have assumed that o, > o,, with strict inequality in general. But 
ifo, =o, =o, the distribution is circular normal, and if we divide all coordinates by 
a, we have a standard circular normal distribution, i.e., with o = 1. We would write 


h * h/o, k= k/o or replace h by oh, and similarly replace k by ok. Also, we would 


divide, u, v, a, b, and R by o. The angle ¢ would be unchanged by the transformation. 


If o, = 9, =o = 1, the case is simplified and the original symbols, h, k, u, v, a, b, 
and R can be retained. 


We would then evaluate a double integral as in (2-12) or (2-13), but call the 
subroutine CIRCV, computing the circular coverage function P(R, d) instead of sub- 
routine PKILL, which computes the elliptical coverage function P(h, k, a,» 9,, R). 
The former subroutine is at least ten times as fast as the latter subroutine. 


This substitution of subroutine CIRCV for PKILL, and the reduction dis- 
cussed below of the double integral to a single integral, also using CIRCV, can be 
done if and only if co, =o,. The reason is that if o, #0, a transformation must be 
made, in order to produce a circular normal distribution, in which unequal factors 
of stretching or compression along the h- and k-axes are used. This would convert 
the lethal circle of radius R into a lethal ellipse and this is contrary to the assumption 
that the bomb’s lethal area is a circle. 


If o, = 0, = 0, involving an initial similarity transformation if o * 1, as dis- 
cussed above, we can even express P, in terms of a single rather than a double 
integral, but at the cost of solving a fourth degree equation in general at every inte- 
gration step to find the points of intersection, possibly four in number, of the circle 
h? + k? = 1? (where r is the variable of integration) with the ellipse. 


Figure Al shows a case in which a circle with center at C and radius r (also a 
circle of radius r + dr) intersects the ellipse in four points. The arcs of the first circle 
(radius r) subtend central angles 6, and 6,. Let 0(r) = 6, +98,. Then the sum of 
the two shaded areas is (6, +@,) 48 dr or 6(r) r dr neglecting infinitesimals of higher 
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Figure Al. Case o, = ,, Ellipse Intersected at Four Points 


order than the first. This area can be described as the area occupied by points that 
are in tie interior of the ellipse and at a distance between r and r + dr from the center 
C of the distribution. Because of the uniform distribution in the ellipse of possible 
positions of the point target, the probability thet the target lies in the shaded area is 
6(r) r dr/(zab), i.e., the ratio of the shaded area to the total area of the ellipse. All of 
these points in the shaded area are at distance r from C, neglecting infinitesinals of 
higher order. Hence, given that the target lies somewhere in the shaded area, the 
conditional probability of its destruction is P(R, r), the circular coverage function as 
computed by subroutine CIRCV. Hence, the total probability that the target lies in 
the shaded area (for a given r) and is destroyed is 6(r) r dr P(R, r)/(rab), and the overall 
kill probability P, is computed by integrating this expression from r, tor,, which are 
respectively the minimum and maximum of values of r such that a circle of radius r, 
center at C, contains points that are on or inside the ellipse. If the point C is on 
or inside the ellipse, r, > r, + U, but if C is exterior to the ellipse, r, > Tt, > 0. Thus, 
we evaluate the single integral: 


1 ft, 
P= — [ 20(r) PR, 1) +r dr (Al) 


tabJT, 
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The determination of tt, and @(r) in general involves the numerical solution of 
fourth degree algebraic equations, but is simpler in special cases as discussed below. 
If C is in the interior of the ellipse, so that r, = 0, we take 6(r) = 2n for all values of r 
such that the circle of radius r (center at C) lies entirely in the interior of the ellipse. 
If C is on the ellipse, we also have r, = 0, but for all circles of radius r > 0, 6(r) < Qn. 


A special case which was programmed and used in checking out the general pro- 
gram of this report is shown in Figure A2. In this case, the limitsr, and r, of (Al) 
are obvious by inspection, and the computation of the function 6(r) of (A1) involves 
only square roots and not the solution of more general fourth degree equations. 


Figure A2. Special Case Programmed 


Here, we have o, = 0, = 1,u = 0, ¢ = 90°, v > a> b, so that the point C is 
exterior to the ellipse. From the figure we clearly haver, =v —a,r, =v+a. It is 
intuitively clear that for v — a <r < v + a, the circle of radius r, center at C, will 
intersect the ellipse in two, and only two points, with equal ordinates by symmetry, 
and this is confirmed by analysis. The common ordinate y of these two points of 
intersection is given by: 


Va? —b*?) (F — b?) + b?v? _ bh 
yo & aa ae a (A2) 
and, when y has been computed, the abscissas of the two points of intersection (circle 
and ellipse) are + Vr? — y?. Hence, the angle 6 (Figure A2) is arc tan (yr? — y?/y), 
and the entire central angle corresponding to the shaded area is twice this or 2 arc 
tan (Vr? — y?/y). Therefore, using (Al), the kill probability P, is given by the 
integral: 


2 vta / 2 yt 
P= — arc tan Le F(R, r) r dr (A3) 
nabJv—a y 


with y at cach integration step, or each value of r, given by (A2). 


Although a program based on Figure A2 and (A2) and (A3), with a> b, was used 
in checking out the general method of this report, we would have a still simpler special 
case if the ellipse were a circle of radius a, with a< v (see Figure A2 and imagine that 
a=b). In this case, the ordinate y of the points of intersection of tiie two circles, 
h? +k? =r? and. h? + (k — v)? = a?, is given by: 

r? + y? — g? 


= Ad) 
y = ( 


and P, would be given by (A3), with y computed by (A4) instead of (A2). 


(b) If a = b, implying that c = | and the ellipse is a circle (in the general case 
where on and o, are not necesserily equal), the positions of the axes of the ellipse 
are indeterminate, and for convenience we can take 9 = 0 so that the x- and y-axes 
are parallel to the h- and K-axes respectively. This would simplify (2-3) which would 
become 


Saat 
k@vtry 


Also, c would be 1 in (2-5). If 0, * o,, we would still have to compute a double 
integral as in (2-12), using subroutine PKILL, but with the simplifications which have 
been noted in the coordinate transformations, (2-3) and (2-11). 


: 
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ALTERNATIVE METHOD OF COMPUTING P, 


We outline here an alternative method of computing the kill probability P,, 
this time using true polar coordinates r, 9 in the ellipse. The curvilinear coordinates 
introduced in (2-5) and also denoted as r, 6 are not polar coordinates in the ellipse 
ifa>b. 


Figure B1 shows a quadrant of the ellipse (see also Figure 1 in Section 2). 
First, we wish to find the coordinate r for a point such as E on the ellipse, with 
given polar angle 6, to establish limits on r in terms of 6. The equations of the line 
OE, and of the ellipse are: 


(B1) 


(B2) 


Figure B1. Polar Coordinates in the Ellipse 


Solving these equations simultaneously, we find that the rectangular coordinates 
(x, y) of the point E are: 


x = (ab/S) cos 6, y = (ab/S) sin 6 (B3) 


r= Vx? +y? = ab/S (B4) 


and, therefore 


where 
: S = Ya? sin?6 + b? cos?@ (BS) 
Equation (B4) provides the upper limit on the variable r in the subsequent double 
integration. 
| 


Now, considering a polar element of area r dr dé (Figure B1) at a point (r, 0) 
inside the ellipse, we must, as before convert its rectangular coordinates (r cos 6, 
r sin 6) to the Chk system, by (2-3). The result is: 


zut = sut + 
h=utrcos6 cos¢ —rsin 6 sing =u +r cos? " (B6) 


k=vtrcosé@ sing +r sin é@ cos¢ =v +r sin(o + ¢) 


By reasoning similar to that used in setting up (2-12), we obtain the following double 
integral in terms of the present true polar coordinates r, 6 


1 2n pab/S 
Px sory | Plu +r cos(é +9), v+rsin(6 +¢),o,,0,,R) rdrd@ (B7) 


The P functions, as in the case of (2-12), are computed by subroutine PKILL. 
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PROBABLE ERROR IN THE WEIDMAN-BRUNNER METHOD 


In the WB (Weidman-Brunner) method as described in [1] and in the present 
technical report, the standard statistical method for determining the probable error 
and confidence intervals is one using the Student t distribution, [4, p. 224], [6, 
p. 178], and [10, p. 189]. For a given problem, c,, o,, R, u, v, a, b, and @ speci- 
fied, there exists a truc kill probability, P, or p,, which we are estimating or approxi- 
mating by the sample mean, p, of a sample of size n (the number of simulated bombs). 
For a fixed and sufficiently large n, these sample means {p} are from an approximately 
normal population with mean P, or p, and standard deviation o (for individual prob- 
abilities p,), the standard deviation of the sample mean B being o/,/n. The point 
which makes the t distribution the appropriate one to use in computing confidence 
intervals is that o (as well as P, or p,) is unknown. The method cited from [4], 
{6], [10], and given in many other statistical texts, is designed for this situation, 
finding confidence intervals for the mean of a normal population with unknown 
standard deviation oc. 


However, in the present technical report, we do not develop the above-mentioned 
method in detail, since it is not used in [1] or in the Weidman-Brunner computer 
program. Instead, we follow this computer program and [1] more closely, and, with 
each sample, estimate the population standard deviation o by means of the sample 
standard deviation. It is believed that this method if followed ur correctly will give 
results that are roughly equivalent to those given by the method which uses the t 
distribution. 


The method used in [1] for computing the probable error (P.E.) is open to 
criticism in several respects. A study of the FORTRAN program given in this reference 
shows that the P.E. is given by the formula: 


n 1/2 
P.E. = 0.6745 | 5° (p, — 5)? al (C1) 


i=1 


where n is the number of simulated bombs, p, is the i-th conditional probability 
determined as explained in Section 1, and § is the mean of the n p,’s, i.e., B is the 
approximation to the kill probability P,. The integer n is usually taken as 25 in the 
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WB method, although the values n = 50 and n = 100 were used in some of the cases dis- 
cussed in Section 3 of the present report. The radical in (C1) represents the standard 
deviation o of the sample of size n, and the constant 0.6745 is a familiar one in 
statistics, where P.E. = 0.67450, since, for a normal distribution with mean 0 and 
standard deviation o, we have: 


0.67450 1 x? 


| 
--— = 1/2 Cc 
ov Qn Fovsrase®? Jax / c? 


20? 
or, half of the items lie between the limits + 0.6745c. 


The denominator under the radical in (C1) should be n — 1 instead of n in order 
to get an unbiased estimate of the standard deviation (6, p. 154]; actually, the use of 
n — | gives an unbiased estimate of the variance o?, and a very small bias in the value 
of o itself. For large n, there is no significant difference between using n and n — 1 
in (C1). 


A second and more serious criticism is that even if we substitute n — | for n in 
(C1), the result is the P.E. for an individual item p, and not for the sample mean p. 
To compute the P.E. for the latter, an extra factor of vn is needed in the denomina- 
tor of (C1) as shown by the following discussion. 


The conditional probabilities p,, p,, .... p, in the WB method are independent 
random variables with identical (though in general unknown) distributions. The dis- 
tributions are identical since, for each of the n simulated bomb drops, the procedure 
is exactiy the same. But we cannot assume that this common distribution is normal, 
since it involves a number of factors including the uniform distribution of the point 
target in the ellipse. 


In this situation, we can get a good approximation for large n, for the variance 
(hence, for the probable error P.E.) of the sample mean p by means of the central 
limit theorem [4, p. 136] and [6, p. 191]. Reference [6] states (changing the nota- 
tion slightly to agree with our present notation) that if p,,p,,..., p, are independent, 
identically distributed random variables, each with mean p, and variance o?, and if 
y, ™P, +P, ++. P, = np, and 
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; then Z, is “asymptotically normal” with mean 0 and variance 1, in the sense that the . 
distribution function F,(p) of Z, satisfies the relation 
if 1 
: lim F (p)= f exp (—u?/2) du C4 
Prats Jin » Pp ) (C4) 


Here, p, is the true kill probability for a given problem of which f is an estimate on 
the basis of one sample of n bombs. We think of a large number of samples of size n. 
The random variable p, takes a large number of values in these successive samples, 
5 with mean p, and variance o* and similarly with p,, p,, etc. 


Since Vn (B — p,)/o has approximately mean O and variance 1, as stated in the 
" formulation of the central limit theorem in [6, p. 191], from which (C3) and (C4) are 
bi taken, B — p, must have mean 0 and variance o7/n and, therefore, B, the sample 
$ mean, must have, asymptotically, mean p, (the true kill probability) and variance 
o?/n, establishing the fact that the variance of the mean is 1/n times the variance of 
an individual p,. 


The facts that the mean of all these sample means is p,, and the variance is 
o?/n, could have been proved by simpler theorems, [6, pp. 135 and 137] but to 
prove that the sample means p are asymptotically normally distributed (and this fact 
is needed in computitg the probable error) requires the central limit theorem. Also, 
if the individual random variables p, were known to be normally distributed, we could 
get very simple and exact, not merely asymptotic, results for the distribution of the 
sample means f. But under the present circumstances, with the p,’s not normally 
distributed, we apparently cannot do better than the present asymptotic analysis. 


rf Since the variance of the sample means p is asymptotically o?/n, as has been 
i shown, the corresponding standard deviation is o/,\/n, where o is the standard 
deviation for an individual value p,. For this 0 we have to use the estimate 
i ([z(p, — BP? /(n - 1)) 1/2, square root of an unbiased estimate of the variance as 
stated above. We multiply this by yn to get an estimate for the standard devia- 
tion of the sample mean f, and by 0.6745 to get the probable error, P.E. (see (C2)). 
Hence, we have - 
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P.E. = 0.6745 [£(p, — B)?/n (n — 1)} 5? (CS) 


Let the P.E. of [1] be denoted by P.E.,, (see (C1)). Comparing (C1) and (C5), we 
see that P.E/P.E.,, = iy n — 1, or 


pe, « Exe (C5) 


P.E.,, is the value printed out for the probable error if the FORTRAN program of 
{1] is run on a computer. These printed values were corrected as indicated in (C6) 
in interpreting the results of computed cases. 


Although the distribution of the sample mean B in the WB method approaches 
normality as the sample size n becomes infinite, there appears to be no simple means 
of determining how fast its approach to normality is, or how large the sample size 
should be to obtain a meaningful estimate of the probable error in the computed 
approximation § to the true kill probability p,. In [9, p. 111), in a discussion of 
this application of the central limit theorem, the following statements are made: 
“In practical terms, a normality assumption for the sampling distribution of X 
becomes reasonable in many cases for N > 4 and quite accurate in most cases for 
N> 10. Hence, for reasonably large sample sizes (4.34) applies to the sampling 
distribution of X computed for any random variable x, regardless of its probability 
distribution function.” Equation (4.34) [9] is 


- N 
(=a) YN =Z (C7) 
Q 


x 


which is stated to have a standardized normal distribution if x is normally distributed. 
The normality is only asymptotic if x has some other distribution. The expression 
on the left-hand side is equivalent to the expression Jn (p — Po Yo of (C3), with 
appropriate changes in notation. Thus, the authors of [9} regard samples of more 
than 10 items as adequately large in order for the use of the central limit theorem 
to give meaningful and even accurate results in most cases. The regular sample size 
in the WB method is n = 25 as stated earlier in this Appendix, and this sample size, 
therefore, appears to be adequate. 


The Chebyshev inequality rigorously assures us, as shown below, that the 50% 
confidence interval determining the probable error will never be larger than the 
interval (B - /20,, 5+ /2 Op), where og (= o//n) is the standard deviation of the 
mean p and o is the standard deviation of an individual item p,. This is true, however 
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small the sample size n may be. This implies an upper bound of ¥2 Oy OF V20/JVn 
for the probable error, which of course, is substantially larger than the asymptotic 


value 0.6745 o//n, for large n, implied by (C5), and based on the use of the central 
limit theorem. ore 


The result stated is derived from the Chebyshev inequality in the following way. 
This inequality, [4, p. 135, (10)], with the notation changed to be compatible with 
the present situation, may be stated as follows: 


vw a5? 
Prob ((B ~ pg! > b)< <= <E (C8) 


where b is an arbitrary real number. Now, we ask what value b must have in order 
that the stated probability shall be 1/2, in order to find 50% confidence limits. Hence, 


we put 
og 
rr) = 1/2 (C9) 
of which the solution is: 
b= V20,= V20//n (C10) 


so that the 50% confidence interval, as stated, is shorter than the interval (5B — v2 os, 
pt ¥2 Ox). Although this result is perfectly rigorous and independent of the sample 
size n, it provides only an upper bound for the length of the 50% confidence interval 
or for the magnitude of P.E. 


For a normal population, the probable error (P.E.), is simply the half-length of 
the 50% confidence interval. More precisely, suppose the population mean is X, q 
in the notation of [10], and the standard deviation (for individual items) is o. Let 
a random sample of size n be drawn, with sample mean X,. Then, the 50% confi- 
| dence interval for that sample is (C,,C,), where C, =X, — 0.67450//n, C, =X, + 
q 0.67450//n, or C, = %, — P.E., C, = x, + PE, since the probable error for the 
4 sample mean in this case is 0.67450/\/n. Thus, the length of the 50% confidence 
interval is C, — C, or 2 P.E., or, the P.E. is the half-length of the confidence interval. 
See the derivation of (C5) above, where the sample standard deviation [Z(p, — B)?/ 
(n — 1)}'/? (for individual items) is multiplied by /Yn to get the standard devia- 
tion for the mean pf, and then by 0.6745 to get the probable error, P.E., for the 
sample mean. 
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We show here how this relation between the P.E. and the 50% confidence inter- 
val follows from [10, pp. 188-189], and we refer the interested reader to this stand: 
ard statistical handbook for further details. 


Let the confidence coefficient be 6, in the notation of [10], and here we are 
interested in the value 6 = 1/2. By interpolating in Table 1X in [10], we find that 
A, OF A, ,, = 0.6745. This well-known constant is the value A, such that 


1 Od, x? 
oVian Sines {[- 21 dx = 1/2 


(se (C2). Then (10, p. 189], C, =X, —A,o//n= 8, - 0.67450/ Vn, C, =X, + 
dyo//n = 8, + 0.67450//n, and the 50% confidence interval is [C, , C,]. The signi- 
ficance of the term “confidence interval” itself is brought out in the following state- 
ment from [10]: “If it is asserted, whenever a sample is drawn, that the unknown 
population mean x lies between the limits C, and C, calculated for that particular 
sample, then the probability that the assertion is correct is 6”, with 6 = 1/2 here. 


In the application made in this report to the WB method, the unknown popu- 
lation mean X corresponds to p, or the true but not exactly known P, , the sample 
mean X, is § or the average of the n conditional probabilities p, for that sample and 
the population standard deviation o in each case is estimated or approximated by 
the sample standard deviation [Z(p, — B)*/(n — 1)]'/?. 


Also, the results quoted from [10] assume a normal population. Thc means 5 
here are not exactly normally distributed, but they are asymptotically normal for 
large sample size n as we have shown, and the normality assumption can be assumed 
to be good for n > 10. 
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APPENDIX D 


DETERMINATION OF POINTS (h,,, k,,,,), (h,,,, > k 


max ) 


OF SECTION 3 


t 


DETERMINATION OF POINTS (hk, .)s as. Kas) 


OF SECTION 3 , 


In Section 3, pages 16-18, we used points (h_,_, K ain» (Hivax? max)» Which were 
useful in simple tests, using very little computer time, for detecting cases in which the 
kill probability, P, is approximately 0 or 1, within a prescribed tolerance. For in- 
stance, if the general tolerance is 0.005, so that only two decimal digits can be assumed 
| to be correct in the values of P, , then, for any case in which Pp. > 0.995, the true 
value is closer to 1.00 than to 0.99. 


In this appendix, we give without proofs, the equations for the coordinates of 
these points Cain? Kmin? Ming x? K nex) (see page 52). 

Suppose we are given (see Figure D1) an ellipse such as those considered in this 
report, with center at (u,v) relative tu tne Chk coordinate system, semi-major and 
semi-minor axes a and b, respectively and major axis making an angle ¢ with the h-axis, 
with -90° < ¢ < 90°. 


k 


(yor 2° Kuor 2) 


(h 


max’? Keax) 


(h. k ) 


vert 2° “vert 2 


h, k 


vert 1° vert 14 


re 0 mr ee en eo re ee ee ee ee 


(Nyor 1’ Kor i } 


min? kin) 


for -90° <  < 90° 


k h 


Figure D1. Values of Hain’ Bmin? Umax? K ax 
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There are obviously two horizontal tangents and two vertical tangents to the 
ellipse. We denote the points of tangency of the horizontal tangents, as in Figure D1, 
as (yor ” Kor 1) and (Dyor 2’ Knor 2) with Knor 2 > Kyor 1 always, i.c., (hy oe 2° 
Kor 2) is the point of tangency of the upper horizontal tangent. The relative magni- 
tudes of h,,, , and h,,. , depends on whether ¢ is a first-quadrant or fourth-quadrant 


angle. 

The points of tangency of the vertical tangents are denoted as (a 1? Kyser Ds 
tn 2’ Kort 2) with Dect 2 > Dvert 1? always, ie., (Hyon 2° Kort 2) is the point of 
tangency of the rightmost vertical tangent (Figure D1). 


tt can be shown by analytic geometry that the coordinates of these four points 
of tangency are given by the following equations. We do not give the proof here. 


(a? — b?) sin ¢ cos ¢ 


Zz u — 
eRe arm nrg ter : 
kor 1 ZY — va? sin?¢ + %? cos? (D1b) 


(a? — b?) sin ¢ cos 9 


ut 
Paor 2 va? sin?¢ +b? cos*¢ 
Kor 2 “V+ va? sin?¢ + b? cos’¢ > 0 (Did) 
Nee 1 ZU va? cos?@ + b? sin? (Dle) 


(a? — b?) sin ¢ cos ¢ 


va? cos’¢ + b? sin?¢ 
bye 2 2¥ + Va? cos*o + b? sin? > 0 (Dig) 


2 _ 2 
vast b?) sin ¢ cos ¢ (Dih) 


k er 
sit Va? cos?¢ + b? sin?¢ 


(Dia) 


(Dic) 


(Dif) 


my ~— 
vert 1 v 
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It is seen from these equations that k,,, , and h. , are necessarily positive, since u 
and v are assumed non-negative and the radicals +a? sin?@ +b’ cos'@ and 


+ Ja? cos?¢ + b? sin? denote the positive root in each case. 


We consider further the particular case shown in Figure D1, in which the ellipse 
lies entirely in the first quadrant. The situation is not so simple if the ellipse extends 
into other quadrants. In Figure Dl, we take h_. #h k.. =k The point 


min vert 1° “min hor 1° 


(Chiin? Kmin) 80 defined has the following property. If (h, k) is any point on or inside 
the ellipse, then: 
h2 Dein (D2a) 
k Pa Kein (D2b) 


Hence, because the elliptical coverage function P(h, k, On» %,» R) decreases with in- 
creasing |hi, and also with increasing |ki, we must have 


P(h, k, Gye My» R)< Phin? ko Oye As R) (D3) 


for any point (h, k) on or inside the ellipse, with strict inequality in (D3), since we 
must have strict inequality in (D2a) or (D2b) or both (the point (h ) itself not 
being on or inside the ellipse). 


min? Kein 


Since (D3) holds for every point on or inside the ellipse, it must hold for the 
average value P of P(h, k, o,, 0,, R) throughout the ellipse. But the kill probability, 
P, is equal to the average value, BP, by (2-17). Hence, we must have 


Px = B = PCa? Kein? On» Oy» R) (D4) 
But the calculation of P, by (2-12) or (2-13) in the general case, for the input values 
used in Table E3 in the present report, requires something of the order of 0.4 second 
or 400 milliseconds, whereas one call to the subroutine PKILL, for computing 
Phin? Kmin? Sn? %x» R), requires only about 5 or 10 milliseconds for three-digit 
or six-digit accuracy, respectively, on the CDC 6700 machine. 

Similarly, in Figure D2, we take hh... hy. ..k.., ™ Kor 2 and by an argument 
similar to the one given above for Chain’ Kin)» We conclude that: 


Py = P> Phi.» Kmax? Sn» %» R) (D5) 


ce ee ce ee Re 2 


The calculation of P(h 
routine. 


kK nax’ [n? %%» R) requires only one call to the PKILL sub- 


The complications which are introduced when the ellipse extends into other 
quadrants than the first are sufficiently illustrated by the case shown in Figure D2. 


The relevant analogs of (D2a and b) here are 


ihi > [h,,,| (D6a) 
ikl 2 |Kental (D6b) 


for arbitrary points (h, k) on or inside the ellipse, and it is seen from the figure that 
(D6a) is not satisfied for all points (h, k) of the ellipse if we takeh,,, =h,,, , a8 be- 
fore. We must in this case take h,,, = 0, since in this case we clearly have h = 0 for 


infinitely many points inside the ellipse. 


a 1? Krert 1) 


Figure D2. Ellipse in First and Second Quadrants 
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The general equaticss for Lr Kain? | k,.,,» Tegardiess of the quadrants into 


which the ellipse exte isis, | 
Semin = max(h,.. ,» 0) (D7a) 7 
koin = max(k,.. ,» 0) (D7b) ‘: 
ha, ™h, 2 >0 (D7c) 
Kas “noe 2 > 0 (D74) : 


with h,,, ,» etc., in (D7a, b, c and d) given by (Die, b, gand d), respectively. Equa- 
tions (D7a and b) refer to the algebraic maxima; thus, for example, max (—3,0) = 0. 


Equations (D7a, b, c and d) are the general equations for all values of the angle 
¢. Only four values of 9, namely, 0°, 45°, 90°, and —45°, are used in the table given 
k. for 


in this report. We give in tabular form below, the values of h k h 


min? “min? 
these values of ¢: 


max(u—a,0) max(u—p 0) max(u-—b,0) max(u—p,0) 


; Kain max(v—-b,0) max(v—p ,O) max(v—a,0) max(v—p ,0) 
h,, ae uta utp utb utp 
k vtb vtp vta v+p 


Here, o = ¥ (a? + b?)/2. It will be noted that the entries here for ¢ = 45° and 
@ = —45° are the same, as they should be. From (Die, b, g and d), it is seen that 
Bert 1 Knor 1° Yvert 2° Knor 22 Which occur in (D7a, b, c and d) are the same for these 
two values of ¢. 
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APPENDIX E 
IDENTIFICATIONS FOR TABLE E3 
TABLE E3 


IDENTIFICATIONS FOR TABLE E3 


Weapon falls in uncor- 
related elliptical normal 
distribution with mean 
point of distribution at C, 
the origin of Chk rec- 
tangular coordinate sys 
tem. Lethal radius of wea- 
pon is R, with cookie- 
cutter lethality function. 
Point target is at random 
point from uniform dis- 
tribution in ellipse. Cen- 
ter of ellipse is at (u, v) in 
Chk system. Semi-major 
and semi-minor axes are 
a and b respectively. Ma- Figure El. Input Parameters o,,0,, U,V, a,b, ¢,R 
jor axis makes angle ¢ 
with h-axis, positive coun- 
terclockwise, ~90° < 9 < 90°, a, = 1 throughout the Table E3; 0 < o, < o,; 
a>b>O,u>0,v>0,;R>0. This is a direct table giving kill probability P, as 
a function of o, (= 1), o,, R, u, v, a, b, ¢. 


Ranges for the variables are as given in Table E1 (also given on page 12). o, 
always has the value 1, a takes the values 0.1, 0.5, 1, 2, and ratio b/a takes the 
values 0.2, 0.5, 0.8, and 1 on every page. Angle @ takes the values 0°, 45°, 90°, 
—45° on every page. Therefore, the only variables whose values are given in Table El 
are R, o,, u, and v. The total number of cases in the table is 59,904, on pages 59- 
109. Three decimal digits are given for the values of P, , but since the tolerance e, 
as discussed in Section 3, was generally 0.005, the third digit may be in error by 
several units, and the values of P, are reliable to only two decimal digits. 
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Table E1. Values of R, o, , u, v for Main Table (E3) 


Pages Values of v 
59-62 0.2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 
63-66 0.2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 
67-70 0.2,0.5,0.8,1 0,0.5,1,1.5,2 0,0.25,0.50,1 
71-78 0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
79-86 : 0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
87-94 0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
(Sse Note) 
95-98 ‘ 0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
99-106] 3 0.2,0.5,0.8,1 0,0.5,1,1.5,2,2.5,3,4,5 0,0.25,0.50,1 
107- 4 1 0,0.5,1,2,2.5,3,3.5,4, 0,0.25,0.50,1 


109 4.5,5,5.5,6 


Note: For R= 2,0, = 0.2: v= 0, 0.25, 0.50, 1.50; 
But for R = 2,0, = 0.5, 0.8, 1: v= 0, 0.25, 0.50, 1 


Examples 


We give two examples of the use of the table. In Example 1 there is only a 
trivial scaling problem. In Example 2 we must interchange the roles of h and k and 
make corresponding changes in u and v, and interpolate with respect to v. See also 
several examples in Section 3, pages 20-21. 


Example 1. A bomber flies east along the h-axis (see Figure El) and aims a 
bomb at the origin C. The bomb falls in an uncorrelated bivariate normal distribution 
with dispersions o, = 250 ft in range, o, * 125 ft in deflection. The lethal radius of 
the bomb is 500 feet. A point target is at a random point from a uniform distribution 
in the ellipse shown, with center at u = 375 ft, v = 250 ft, semi-axes, a = 125 ft, 
b = 62.5 ft, with the major axis inclined at 45° as shown. Find the single shot kill 
probability P, . 
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Solution. We must normalize by dividing all distances by o, 7250. The unnor- 
malized and normalized values of the input variables are as follows: 


On page 89 of the table, entering the table with the normalized values of the variables, 
we find that P, = 0.540. Since the tolerance ¢ is 0.005, the last digit is in doubt, but 
the figures 0.54 can be regarded as reliable. 


Example 2. Suppose we are given the values o, = 160, o, = 200, etc., as on the 
first line of Table E2. 


Table E2. Data for Example 2 


‘ Original Data 
2. Normalized 


(divide by 200) 

3. Rearranged with 
O. > o, 

4. Tabulated Value . . B 0° | 0.234 
(page 84) 


Tabulated Value 


By interpolation, P, = 0.200 for data on line 3. 


58 


fe flat Pe eee te we pe Ne eS OM Se UP Re Rey ee Cored oe EOS Ske Pe ae ee eee eg Re ee ee ot dea ny a ee ee TN ae ee ee te a ee Re tee ee eg ee! ee ee Ore Ue ae ee ee 


We first normalize by dividing distances (all variables except ¢) by the larger 
standard deviation, o, = 200, producing the data on line 2. Then we must rearrange 
with o, 2 o,, in order to use the table in this report, producing the values on line 3, 
The values of o, and o, are interchanged, and, as a consequence, also the values of 
u and v are interchanged from the values on line 2. 


(a) (b) 


Figure E2. Normalized Data for Example 2 
(a) Line 3, (b) Line 3 of Table E2 


However, a and b are still the semi-major and semi-ininor axes respectively, of the 
ellipse, so are not interchanged, but a consequence of the rearrangement is that is 
changed from 90° to 0°. A consideration of (a) and (b) of Figure E2 will convince 
the reader that these figures represent the data on lines 2 and 3 respectively, and that 
the cases are equivalent. 
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However, the value v = 0.7 is not given in the table and so we must interpolate 
between the values v = 0.5 and v = 1 which are given (see lines 4 and 5). The results 
from the table, page 84, are P, = 0.234 for v = 0.5 (line 4) and P, = 0.149 for the 
data on line 5, with v = 1. Interpolating linearly, we find that for the data on line 3, 
P, = 0.200. Here the linear interpolation over a relatively large span, v = 0.5 to 
v = 1, introduces an additional error which is difficult to assess on theoretical grounds. 
However, this case was computed directly by the DJH program, using the original 
data (line 1, Table E2). In the computer program, o, is not required to have the 
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value 1 or to be not less than c, , although such restrictions were imposed in computing 
Table E3 for this report. The direct computation, with a tolerance ¢ = 0.005, gave a 
result of P, = 0.210897 --, rounding to 0.211 and thus differing by slightly more 
than one unit in the second decimal digit from the result obtained by interpolating in : 
, Table E3. Since the overall tolerance was ¢ = 0.005, the interpolation introduced 
only a small additional error in this case. 
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APPENDIX F 


FORTRAN LISTING OF DJH PROGRAM 


Preceding page bank =? 


ST een ee we 


ee ee 


4 a an ate ae 
eee eee eee - - 


ee LNT A gn eae = 


FORTRAN LISTING OF DJH PROGRAM 


EXPLANATION OF FORTRAN LISTING FOR P, 


Routine LDBSMP evaluates the integral 


nf" f P(h, k, o,,0,, R) r dr dé 
with 
+R pht VR? -(x-h?P 2 2 
Pvh, k, Oy: Oy» R) “f{" tS exp _! Lome Ae dy dx 
k— JR? -(x—h)? 2 \o,? 0? 


evaluated by subroutine PKILL, where 


h=u+rcos 6 cos¢ — cr sin @ sin @ 
k =v +rcos0@ sing + crsin 9 cos@ (See (2-11)) 


Identifying the other variables 


is major axis of the ellipse 

is minor axis of the ellips: 

is b/a 

is the abscissa of the center of the ellipse along h axis 

is the ordinate of the center of the ellipse along k axis 

is the standard deviation in the h direction 

is the standard deviation in the k direction 

is the lethal radius of the weapon 

is the inclination of the major axis of the ellipse to the h axis measured 
positive counterclockwise from h-axis, —90 < ¢ < 90. 


ES 28: SE Or a 


Routine LDBSMP calls subroutine ARIEL in the followin3 manner 


CALL ARIEL (u, v, a, b, o,,0,,¢, R, P, , tolerance) 


with tolerance set to 0.005. Its basic function is to set up the input and locate the 
output. 


LT 


Subroutine ARIEL calls the double integration routine DBSMP to compute P, . 
When o, ™ 9, the routine is much faster because the circular coverage subroutine 
CIRCV can be used instead of the elliptic coverage function subroutine PKILL, (see 
page 32). 


Subroutine DBSMP is called by 


CALL DBSMP (xu, xf, yu, yf, eps, N, g, PBXY), 


xu is the upper limit of integration of x 
xf is the lower limit of integration of x 
yu is the upper limit of integration of y 
yf is the lower limit of integration of y 
eps is a tolerance which is used to terminate the integration, when 


[,, —4a,| << epal(l + yu ~ yf) (See page 23.) 


where 


|, approximates P, with partitionings of the (x, y) (or r, 6) intervals 
of integration and n, is a refined partition of n, . 


N is the maximum number of subdivisions of an interval of integration 
in either the x or y direction. 


g is the location of the value of P, . 


PBXY is a subroutine which evaluates the integrand at any point in the 
interval of integration. 


Subroutine INLINE is used by D3SMP. 


Subroutine PBXY—3 replaces PBXY in DBSMP when o, =o, ,i-¢., when CIRCV 
is used instead of PKILL. 


Subroutine CIRCV is called in the following way 


CALL CIRCV (R, d, 1, P, IR), 


115 


_ ORRIN ERR eee 


where 


R__ is defined above 
d sis the distance from the origin to a point of integration (h, k), i.e., . 


d= Sh? +k?. 
P is the location where the value of P is stored 
IR isan ertor indicator which is not applicable here. 


CIRCV uses the special funciions ERF and ERFD in tabulated form. They are 
defined as 


ERF(t) = — f ter? ay 
va JO 


ERED = —2- er!” 
nT 
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10 


15 


20 


25 


30 


35 


1101 
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PROGRAM LOSSMHP COUTPUT 6 
COMMON/ONLY2/FR(1025) sEPSLONPTS oN39N5.N7,I0PS5 
TOPS2i 

FORMAT ¢€ 1H0,6€22.15 ) 

FORMAT ¢€ 1H1 Y 

PRINT 97 

PI=3.214615 92653 58979 

UL=16f. 

Viz&60. 

Aiz200. 

83-100. 

SxXi3460. 

Sv3=200. 

P32(90.°PT)/180. 

R3=300. 

E°S12.0002 

CALL ARIEL € U1Le W191 983 eSX1LeSV3eP 3, R3,GS,EPS1) 
PPINV 91,U1,V1,A1,83 4SX1 oSV¥3,P3,R3,6S 

CALL EXIT 

ENN 


SUBROUTINE ARTEL (UUs VV 9 Ae M pSIGXySIGY, ANGLE, ARKILL»GS,EPS1 ) 
EXTERNAL PSXY ,PBXV 

cu 4T101/7U 5 Vp PHI, CEE p SIGMAX y SIGHAY gRKILL, SNPHTCSPHI 
2 

vavv 

RKILLZARKILL 

SIGMA YsSIGV 
SIGMAX=sSIG£ 
PI=3. 146159265 358979 

CEEZSAZA 

TMOPT 22. *PT 

PHIZANGLE 

SNPHTISS IN CPHL) 

CSPHI=COS (PHI? 

xUaTwoPT 

EPLZE PSL PIFAFA 

E1sEP 

XL 20 

YUSA 
Y¥L20 

TF € SIGX.EQ.SIGY) GO TO 1102 
CALL OBSMP (€ XUpXL p¥Uy YLyE 1) NMAX,GS,PBYXY ) 
GSz(2.°GS) S(TWOPIFAFA ) 

GO TO 1154 

LONT T NUE 

UsU/SIGX 

VaV/SIGX 

SIGMAX=1, 

SIGMAY=1. 

RKILLSRKILLZ SIGX 

YUZA/SIGX 

ELZEPSA PIS y UYU 
CALL DBSMP ( XUgXL,VU, VL, EL, NMAX,GS,PBXYS ) 
GS2GS/C PI*i CA*AD/SISIGK*®SIGK) ) 2 
CONTINUE 
RETURN ‘ 
ENO 


we RE fe eet ee vane name oe 


SUBROUTINE OBSMP(XUP,XLO, YUP sYLO,EPSR, NMAX,EVALFN, POXY ) 
CONMON/ONL Y1/FR(1025) ,EPSL»NPTS »N3_N59N7,TOPS 
PTI2#23.1415 9265 3589 79 

NMAX 31025 


| 
| 
| 


EPSL SEPSR/ (1 «+ (YUP<YLO)) 
Hz (XUP-KLO)/2 6 


; 10 c 
c FOR 3 LINES 
Cc : 
1Y=0 
: Ne3 
15 HK2( YUP=YLOD /2. 


00 100 JJei,3 
TYsIVet 
; FY2YLO+FLOAT(E JJ-1) FHK 
100 GALL INLINE( FY, Judy TVs POXY, H » XLOs NHAX ) 
20 CONSZSFR (1) +FR(3) 
OL VAL # (HK* (CONS2+4.°FR(2)))9/3. 


SAVE LINES ALREADY INTEGRATED AND MULTIPLY OY 2 


aad 


25 150 PLOFn, 
DO 200 Je3yN 
LaneJe2 
NUMs 2"Lo4 
FR(NLM) SFR CL) 
30 200 PLO=PLOSFR (L) 
PLO=2.*PLO 


COMPUTE ADOIVIONAL LINES AND MULTIPLY BY & 


35 PLN®0. 
N=2*Ned 


QOD 


; HK=HK/ 2, 
' ODO 300 JJ22,N,2 
TVeTy +1 
| 40 FYSYLO*FLOAT( JJ-1) FHK 
CALL INLINEC FY» Judy ITV_ PBXVY, 4 9 XLOp NMAN D : 
$00 PLNMPLNCFR (JU) 
PLNS4e*PLN 


4&5 INTEGRATE RECTANGLE 


EVALFN#(HK*® (CONS2 ¢PLO*PLN)) 7/30 


TEST FOR DESIRED PRECISION | 
TF (ASS (EVALFN=OLVAL) <CEPSR*AGS (EVALFN))) 500,500,400 
TF CABS (EVALFN=OLVAL) (EPSL ») 500,500,400 
400 OLVAL=EVALFN 
IF ( NeLToNHAX) GO TO 150 
C 500 RETURN 
55 500 CONTINUE | 


IF € I10PS.EQ.1 ) EVALFPNSEVALFNO(61./250) 9CEVALFN©OLVAL ) 
RETURN 
ENO 


Q0o0o ana 


REGISTER ALLOCATION 
& REGISTERS ASSIGNED OVER THE LOOP GEGINNING AT LINE 26 


f 
q 
{ 
{ 
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SUBROUTINE INLINE( YVoJJy IZy FOXY, OHOy XLO» MMAX ) 
COMMON/ONLY1/FR(1025),EPSL»NPTS sN3_N5 »N7y TOPS 


INITIALIZATION FOR FUNCTION EVALUATION AT 3 POINTS 


IvzIZ 

H=OHO 

Nz=3 

IXzi 
; 10 X=XLO 
; . CALL PBXY( X_ Y»p IV, IX, ANS) 
7 CONS12ANS 

1X22 
t X=XLO+H 
| 15 CALL PBXYC Xy Yy IV, IX_ ANS) 
[ 


iv kre ye 


SUM2= ANS 

1X23 

K2=X+¢H 

CALL PBXY( Xp Vy IV, IX, ANS) 
i 20 CONS12=CONS1¢ANS 

OFR#H* (CONSL +4 o*SUM2) 7/36 


SAVE PREVIOSLY COMPUTED VALUES AND MULTIPLY BY 2 
25 25 =PFOR2.*SUM2 


COMPUTE NEW VALUES AND MULTIPLY BY & 


wet Pee ati Fea 


aan aan 


SUMG=0. 
; 30 Mx 2° Mot 
H=H/2 
DO 30 J=2s¥;° 
IKzIX 41 
: X=XLO+FLOAT(CJ-1)*H 
35 CALL PBXY¥( X»y Ys IVs IX ANS) 
t 


a a 


30 SUMG=SUM44ANS 
PFNZT4 o*SUMG 
SUH2 = SUM2 ¢SUMG 


i Sala 


INTEGRATE ALONG LINE AND TEST FOR DESIRED PRECISION 


wo 
aon 


| FRC JJ) =(H9 (CONSS PF NtPFO)) 73. 
IF ¢ MeEQeS5 } GO TO 50 
Cc IF CAISC(FRCJJ) ~OFR) -EPSL*ABSCFR(JJD)0) §=409 40, 50 
45 IF CABS (FR (JJ) ~OFR) -EPSL » 40,40, 50 
50 OFR2FR(JJ) 
IF(CM.GE.MMAX) RETURN 
IF (LY eNE ot eANDeIXSEQeNPTS) L¥s1 
GO TO 25 
50 4&0 IFCIYeEQ.1) NPTS=IX 
N7=N74M 
IF ( N3-LT.M ) N3SM 
NS5S=N541 
IF € IOPS.EQ.1 ) FRI JID SFRO JID 402/12 5,99CF RISD -OFR ) 
55 RETURN 


END 
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wu Fa #N 


SUBROUTINE PKILL (RRR, SAR, SYR, bHR,RKR, RESULT) 
OIMENSTION WN(6) XI (26,6 DoWI(B G6 ) 
THIS PKILL IS FOR 3 OLGFT ACCURACY. 


REAL K 
INTEGER OG 

DATA WNHE4/.99999975/» TONMES/1 69999995 /,A/5, 652 /,ASTAR/S S747, ALS 
Ae 8927 AS /5..233/9SQRT2/ 1.41422 S6/,NN/6y 8912916, 20 26/ sy THSQPI/101283 
8792/7 »SQ202/.70730678/, TOL AE-07/, 

CXI/= 6 93246951 9-2 661209399~ 02 38 61919,.238619195 66120939, .93246951, 
D1870..05— 296028986, — 2796666469 =. S25S532h19 ec 1BTHSEGH, 518303864, . 5255 
E3244 » «796666485 6 96028986, 26°00 09= 2 981560635— 2906117265 —. 76990 267, - 
F 58731795 9= 636763150 y— 0 12523941). 12523 3419 o 36783150 p05 87317955 0769 
CIOZGF » 2904117265 « 98156 063, 1270 op — 2 9894 0093, — . 94H87502,-°, 86963120, 
He 6 75540441 9-0 61787624, ~. 4S 8026 785 > 128160355 » - 095012510 F-01,.950125 
I10E-01 5228160355, .458601678 5 661787626, . TSSHOG G1, . 86563120, . 94457502 
Jee 989400 93,890.09 —o 99312 860) — 0 96397193 9-0 PL ZZIHHS» — 0 83921697, >. 7%6 
K 331915 -. 63605 368, — 651086780 y~. 37370609, —. 227 785855 -. 765265 21 E015. 
£76826 521 E-025 0227785859 «37370609» «51086700, . 63605368,.74633191,.83 
M9116 97 9e 91223 443, 6 963971935. 99322860 ,4% 0.09 - 699518722, —-. 97 472856, 
Ne 93627435 9—. 8 8648 5539208200 0199p >6 740124199 - 64809365, -. 56542167, ~ 
06433793519 —- 231504268) = 2191118867, —-. BH0S669SE- 01, . 64056893 E-01,.1912 
PLBOT 5 o 3150426 8 9043379351 po SUSH2L HT » 66609365 y .74012619,. 8200019956 
28 8641553, .938274655,.974 72556, . 99516722/ 

DATA WI 017132449 5 » 560 7625 790% 6791393, 66791593» 6360 76157, 01713246 
W9918*%0.0pe108 226545 .222381035 «31370665, .36268376,.36268378, «31437 
80665 » 02223610 3, 6101228545169 0. Oy oh 717S33GE“O 1, . 106939335 16907 B33, 
Co 2OT1G 7H 39 FS3H9Z5 Gy 249187055 o2HU914705 5 23349254, . 20326763, .16097 
0833, 610693933, .47175S36E=-81, 129000» -2P7ISZHS9E—01» 622535 PHE~ 01 095 
ELS651 26-015 2124628975 .1495 95995 16915652, LB 2603H2 919945061 ,0189% 
FSOb1» 018 260542 5 0169156525 6 1495 9599 5 022462897 9. 95158S12E 901 » 622535 
G24E-02 9 0 27152 459E $01, RF 0.0, L7OLGDOTE M1, HOCOLG SOF 15. 62672048E=2, 
He 6327674 2E 915 0 LDA 99012» 0 116194535. 1316 BBG, . 142096115. 1691 7299 5.15 
1275339 5015275 339» 0169172999 0142096115 0 131688 Ge yo 11629453 » 220193012 
Jy e B3Z27E7HZE “1, 6267206 BEML». 48 6014 30E—1, 176 108 O07 FHi4 470.0, 123612 
K3OE< 02 pe 28532 SE9E O18 pe HH2T7HS9E= 04 yo 89298585 F-00157 S34GNB1 E045. 86 
£19016 2E-01 5.97618652F-0 25610 74 Hb 27» 082550567 5 0122670475 01258 3746,, 
M12793620 50127938205 oAL2ZSHIT by. 12167047, h1S5S05675 10727, . 976186 
NS2E=0 159286190 2628-019 oe 7SIHGLBLE 01» 5929858SE—01 5 4427 7H39E=-01 5.287 
0531389E-01,.12341230E-01/ 

RaRRR 

SXsSKR 

SY*SYR 

Ha#HHR 

KERKR 
H2aH * HK 

IF (SX-SY) 192,53 
S3SY 

GO TO 4 
H = SQART(H2) 

K=0. 

S=Sx 
TEMP 2R-A5*S 
IF (TEMP) 596,6 
TEMP#ReAL*S 
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ero. 
EF ote 
EF OOS * 
FF 80% 
FFOOS 
EF O06 
EF 007 . 
FF 008 
EFeO9 
EF OL8 
FFOAL 
€F 042 
EF 013 
EFO14 
EF oas 
EF O86 
FFOL? 
EFOL8 
FFO19 
fF O20 
EF@2t 
EF 022 
EF O23 
EF O2s 
EF O25 
EFO2 
EF 027 
fFO2A 
EFO29 
€F 030 
EF O32 
EF O32 
EFOSS 
FFO34 


70 


7s 


90 


ge, ee aE ye ee ee ee ee ee ee ee ee 


100 


105 


2°06 


27 


26 
29 


TF CTEMP* TENP-H*H-K*K? 85859 
IF (CTENP* TENP“H*H@KFK) 557,57 
Psi. 

GO TO 75 

P20. 

GO TO 75 

TF (R-H*OAL*SX) 826510 
TF(R-K+AL*SY) 6585924 
SX225X*SQRT2 
SY22S¥*SQRT? 

T2H* SV 4K *SX-A1 *SX2°SY 
ITFCT) £3513,12 
TF(R*R*(SK*SK +SY*#SYD TT) 8,8, 143 
ASXsA*SX 

ASY=A*SY 

TEMP ZR -H “ASX 

IFCTEMP) 14514,15 

GO=0. 

GO TO 16 

GO=STEMP/R 

RARY2R-K=-ASY 

IFCRARY) 17517,18 
GOP=0. 

GO TO 219 

GOPZR ARY/R 

TF CH-ASX) 20,20,21 
Gizl. 

60 TG 22 

Giz (R-HFASKI/R 
TEMP=K-ASY 

IF(TEMP) 23,523,246 
GiP=1. 

GO TO 25 

GLPs(R-TEMP)/R 
¥I#G61-G0-G61P+60P 

IF CVI) 27,27,26 

RARY=SX 

TEMP2h 

Tz 63 

UsGu 

SxXzSy¥ 

Hak 

G6is61P 

GO=GOP 

SY¥#RARY 

K=TEMP 

GiP=T 

GOP=y 

GiGOP2G1*G0P 
61PG0#G1P*G6O 

TF (G16 OP -G1iPG 0) 32 ,52,28 
TFC ABSC(CVI) = TOL) 314,31,29 
U2604+61+G0P4+G1P 

T2Z.* (GO*G04G 2°64 6 OP* GOP +61P* G1) -UFy 
IF(T) 30,32,32 
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q 
i 
\ 
t 
\ 


£15 


420 


125 


130 


135 


140 


145 


150 


155 


160 


165 


30 
31 


40 


a1 
42 


43 
4&4 


45 
50 
51 
52 
53 
54 


55 
56 


57 


$S 264 .°G1G4P*G1PG0 
TE (SS<T*T) 34, 32,32 
RARY=SX 

TEMP=sH 

T2G1 

UzG0 

SX2SY 
42K 
61361°P 
GO260P 
SYZRARY 

K= TEMP 
GiPsT 

GOPsY 

Ei = SQRT(Gi) 

EO = SART(GI) 
EN2(E1-E0) *(.27"*R/SX#(R/SY)D/1.08"% ABS(R-K) /SY¥#10.)) 
SX2eSX*SQRT2 
SY2eSY¥*SQRT2 

AS YzA*SY 

IF € EN-200 ) 4&4540,40 
OG=4 

GO TO 50 

IF € EN-o7 d 43942942 
0G6=3 

GO TO $0 

IF ( ENo.S ) 455 4b, 4h 
0G=2 

GO TO 50 
0G=4 

TA=1 

Te=41 
IF (CH) 529512,52 

TH=2 
60 TO 53 

TH=1 

IF (K) 56,555,546 

TK=1 
GO TO S56 

IKs2 

P=0, 

N=NN (OG) 

00 74 I219N 
TIJUs CEL“EOD F(X ICI, 06)41.9/2-4+E 0 
V=TIJ*TIJ 
Y=R* (1 = V) 
XL2V¥*(2.-V) 

XQ s SQRT(KL) 
KL=XQ-(ASYV-K) 7R 
Q=XQ= (ASY*#K)/R 

XI Jz (H-U) /SX2 
Xz~-XIJ*XIJ 
Fl = EXP(X)* THSQPI 
GO TO (57,58) ,TH 
VI Js (HOU) 7SX2 
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170 


175 


180 


165 


1906 


195 


59 
60 


61 
62 
63 
64 


65 
66 


67 


68 
69 


70 


71 


72 
73 


7h 
75 


XaoV¥TJ*VIJ 

F2 2 EXPCX)* THSQPI 
FIJsFi+2 

GO TO (59,72,69) ,T@ 
FIJ22.°F 1 

GO TO (59,72,69) ,IA 
IFCXL) 60,71,72 

0112 (K+R*XQ)/SY2 

IF (012-A"%SQ202) 62,61, 61 
OLIZASTAR*SQ202 

Pt = ERFNWL(044) 

GO TO (64,70) »TK 
001=2(K-R*X0)/SV2 

TF (0014A"°SQ202) 65,65.66 
00 12-ASTAR*SQ202 

P2 = ERE NNL (002) 

ITF (O02) 67,67,68 

22 Pi+ ABS (P2) 

GO TO 69 

2=P1-P2 
P=PewWl (I ,0G)*FIU*®Z*TIU 
€Q TO 7% 

P22P1 

60 TO 67 

PLsWNMNEG 

TAa2 

T8622 

TF(Q) 63,73,73 

Z2 TONMES 

TAz3 

Tsz3 

60 TO 69 

CONTINUE 

P2P*¥R* (E1-E0) 7 (SK2%4,) 
RESUL T=P 

RETURN 

ENO 
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SUBROUTINE CIRCV (RK, OCy KR, PV, IR) 
rR = 0 
TF (28K) 109 205 30 
10 IR= 7 
5 15 PY = -4. 
RETURN 
20 Pv = 0. 
RETURN ane 
40 IR § 
[ 40 GO TO 45 
100 PV = i. ae 


RETURN 


C= 0. 


Cc ‘ 
30 A = RK } 
15 8 = ABS(DC) 

IF (KR) 40, 220, 50 


Cc 
S50 WLOC = A-B 

20 IF (WLOC + 661094103) 20, 20, 90 
' 90 IF (WLOG = 6.54466793) 1120, 100, 100 
j 110 3A = A*S 
KFLAS = 0 
X BS SIGN(1e,WLOC) 
2s IF (BA = 100) 220, 120, 240 


Cc 
220 IF (A = 625446793) 2225 100, 100 
222 IF (3) 2905 290, 226 
224 OFF = AYA 
( 30 SA = BR 
IF (COFF/SA = OFF = 40.) 230, 230, 290 
230 OFF = A/8*.5 
A = (1.-B) *OFF 
B= (4.¢8) *OFF 


ne ee, = 2 
ee ee SOE nee nae 


| 35 KFLAG 2-1 

| ga = A*B 

WLOG = A-8 

| C 

| 120 S © ERFD(.7071067 6" A) *ERFD( «70 71067898) 
40 DE = ATA + B86 


DD = (2./DE*5A) 992 
Q = ABSCWLOC* (A+B) 70E) 
20 = .78539616"S 
WO = Q-Q°720 
45 z2nm2 = 20 
W2NM2 = WO 
STG8 = QO. 
IF (KFLAG) 130, 160, 160 
c | 
50 130 GC = Cet. i 
IF (C=1000.) 140, 250, 150 
440 Z2NM2 = {.5/C*B3A)**2 * Z2NKH2 
OFF = O0*W2NMZ 
W2NM2 = COFF = oS/CFOFF) = (ee /DE"C + 10) *Z2NM2°Q 
55 SIGS = W2NM2+SIGB e od 
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70 


7s 


100 


105 


110 


159 
160 
170 
180 


190 
200 


240 


270 
280 


290 


IF (ABSCW2ZNH2) = .0000005) 154, 150, 130 
PV = WO + SIGB 
RETURN 


SIGA = 0. 

C = C+i. 

IF (C=-21000.) 190, 200, 200 

Z2NM2 = (BA/C*.5)9*2 * ZONK? 

SIGA = Z2NM2+SIGA 

OFF = DO*W2NN2 

WZNM2 = (OFF = .S/C°OFF) = (a. /DE®8C © 1.) *272NN24Q 
SIGS = W2NM2 + SIGB 

IF (AGSCW2NN2) - .9000005) 190, 190, 170 
IF (ABS(Z2NM2) = .0000005) 200,200,170 
OFF = 1.-Z0-SIGA 

PV = ((W0 + SIGB)*X + OFF)4.5 

RETUON 


OFF = ABS(.70710678*WLOC) 

SA = ERFOCOFF) 

$B zs -ERF(OFF) + 1. 

S = &e"BA 

SQ = SQRT(2."° BA) 

SIGA = SQ*SA 

W2NM2 2 .25°SA/SQ 

Q = (A+B)/S0" S8*, 70710678 

SI6B = Q 

DE = ASS ((CAtB) /S* WLOC) 

DO = WLOC/S*®WLOC 

C = Cle 

IF (C-1000.) 260, 280, 288 

WLOG 2 Q*00 

Q = DE*W2NM2 + .S/C*WLOC - WLOC 

SIGB = Q # SIGB8 

Z2NN2 = (C42. = 1.) *W2NN2 

SIGA = Z2NM2 + SIGA 

W2NM2 = (2.540) /(4.¢C) *Z2NN2Z/5S 

IF (ABS(Q) = .0000005) 278, 270, 250 
IF (ABSC(CZ2NM2) = .0000005) 280, 289, 250 
Py = ((1e-SISB)"X = SIGA/S + 1.) * 25 
RETURN 


OFF = .70710678°A 

SA = ERFOCOFF) 

SB = 1.-ERFCOFF) 

OFF = 1.-B*8 

WLOC = SQRT(OFF) 

M2NM2 = .35355339/A*B/WLOC 
Q = SB/WLOC 

SIGB = Q 

SA « SA*8/ OFF 

$8 = 8*B/0FF 

OFF & SB/A*#2 

C = C44. 

IF (€=-1000.) 310, 320, 320 


125 


{ 
l 


oot moter ne 


} 310 WLOC = Q*sBe 
; : Q = SA*H2NM2 + MLOC*.5/C = MOC 
SIG8 = Q¢SIGB. 
; W2NN2 = (2.%Co10) @HZ2NM29OFF*(C4.5) / (Cote) 
115 IF (ABS(Q) = .0000605) 320, 320, 300 
320 PV = 1.-SIG8 
RF TURN | 
END 


‘ 
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10 


15 


20 


2s 


30 


40 


45 


50 


55 


FUNCTION ERF(X) 


DIMENSION TERA(425), TERD (e521) 


PODBVOFOWSUN KP OGAWA RNFAaANnNe 


OATA (TERA(T), 


BN OeErNOREUNPH PU BSVNOWEFUNKDOANFYW FAHY 


DATA (CTERACI)» 3245125) / dep 13°. 99999999, 6%. 99999998, 


5* 299999997» 3%. 99999996, 2° .99999995», 2% 9999999, 
2° 69999999Se 2°.9999999Z, » 99999992, 99999990, 2%. 99999989, 


099999988, 
o 99999982, 
0 99999976, 
099999962, 
© 99999905, 
¢ 99999920, 
099999885, 
99999837, 
©99999766, 
099999672, 
099999540, 
299999356, 
¢ 99999203, 
© 99998757, 
299996286, 
09999 7647, 
© 99996787, 


0999947605 
99992931, 
099990508, 
« 99987316, 
©9998 3132, 
099977671, 
©9997 0584, 
«99961429, 
© 999469665, 
©9993 4620, 
099915478, 
© 99691242, 
99860712, 
0 99822446, 
09977H7155 
099715485, 
°99642346» 
099552685, 
© 99442628, 


DATA (TERA(T), 


©99215622, 
0990346861, 
98617420, 
© 98557850, 
©962469279, 
0 97 686284, 
097456701, 
«© 96981622, 
© 963635407, 
e 95665725, 
«9690 1604, 
«94001503, 


2 89999957, 
299999964, 
2 99999972, 
99999959, 
099999940, 
099999914, 
99999877, 
¢ 99999825, 
«99999751, 
© 99999649, 
o 99999507, 
99999311, 
o 99999062, 
© 99998674, 
09999617 3, 
¢ 99997495, 
99996582, 


12126,223) / 199995632, 


0 9999RK SH, 
« 99992699, 
© 99986993 6, 
© 99986567, 
o999862151, 
© 99976397, 
«99966934, 
2 99989305, 
099946942, 
2 99932149, 
099911073, 
2 99885682, 
099853728, 
© 99613715, 
099763661, 
299702083, 
099625810, 
0995322275 
099617933, 


T2Z26,321) 4.9930899G, 199279163, 


©99182075, 
298996516, 
098769094, 
© 98500283, 
9818 LVRS, 
« 97803809, 
097360263, 
¢ 96642350, 
2 96237290, 
o9S537616, 
09%731240, 
¢ 93806526, 


© 99999986. 
2 99999979, 
oe 99999969, 
e 99999956, 
0 99999936, 
oe 99999908, 
¢ 99999668, 
099999812, 
¢ 99999734, 
¢ 99999624, 
099999873, 
«99999264, 
© 99998977, 
o 99998586, 
099998053, 
099997333, 
2 99996368, 


099994090, 
2 999920462, 
© 99989335, 
2999865776, 
099981128, 
0 99975056, 
099967198, 
099987071, 
099944083, 
o 929275 06, 
©9990 6457, 
29987 9861, 
0998646423, 
© 99804894, 
099752529, 
© 9968 8046, 
09960 858i, 
0995211161, 
099392257, 


© 99147249, 
oe 98952454, 
© 98719028, 
« 964460701, 
o 96140692, 
097720684, 
o 97262812, 
©9672 7678, 
«96105351, 
oBSA5244, 
0 I4SS61bG, 
0 93606312, 


299999985, 
99999977, 
o 99999967, 
- 99999982, 
« 99999931, 
« 99999901, 
«99999858, 
99999798, 
» 99999715, 
999995986, 
0999998 36, 
099999214, 
© 9999890 8, 
99998492, 
99997926, 
299997162, 
999961 36, 


099993725, 
99991559, 
099968699, 
099904961, 
99960029, 
099973640, 
0999653741, 
099954724, 
099961060, 
99923666, 
99902620, 
099873766, 
099838763, 
99795065, 
099740702, 
099673461, 
099590635, 
99469205, 
099365565, 


©99112103, 
«9890905 0, 
«98667167, 
+ 98379046, 
-98037559, 
97634838, 
097162273, 
296618515, 
. 95969503, 
-95226512, 
+96376220, 
293400796, 


299999983 » 
2 99999976, 
© 99999964, 


099999988 , 


© 99999926 » 
2 99999893, 
© 99999848, 
09999976b, 
© 999996946, 
© 99999570, 
«99999397, 
» 99999160, 
«99998835, 
299998392, 
099997791, 
099996980, 


099995890 / 
© 99995067, 


© 99993339, 
099991048, 
+ 99988026 
09998060, 
© 99978881, 
099972152, 
099963450, 
299952287, 
© 99937928, 
99919679, 
. 99496551, 
099867387, 
+ 99830795, 
099785114, 
09972361, 
099658215, 
© 99571945, 
099666372, 


099337823 / 
09924 7932, 


© 99073595, 
2 98864055, 
« 98613659, 
998315259, 
o 979621 78, 
© 97546202, 
© 97056569, 
0964 89786, 
© 95829657, 
© 95067330, 
094191372, 
o 93189863, 


60 
65 
70 
7S 
80 
85 


| 
i 
\ 
i : 


95 


100 
105 


ANUS anre POAnTNFCUNKP CBOVFVrunrD OSVNAREFUNYEY BDOVNFVS ANP KP WeFNFus 


© 9297 34195 
«91805010, 
2906863743, 
« 88997067, 
087332616, 
© 85478621, 
283423150, 


DATA (TERACI), 


«79690821, 
077066806, 
74210096, 
0711315563, 
067780119, 
264202933, 
260385609, 
056332337, 
252069968, 
047548172, 
e42839236, 
037936205, 
032662676, 
02763 2639, 
022270259, 
016799897, 


092751363, 
292553388, 
090280040, 
« 66678789, 
086977330, 
«880834602, 
282987029, 


T#322,415) / 281156356, 


079184325, 
«76514271, 
073610345, 
«70467606, 
67084006, 
063458583, 
59593650, 
255493925, 
051166826, 
«46622512, 
» 42873870, 
«36936453, 
¢ 32828350, 
226570006, 
© 22183989, 
015694703, 


«925235994, 
291295551, 
269909620, 
« 88353301, 
« 86614359, 
«64681050, 
0 867542365, 


078668732, 
© 75952376, 
073001043, 
© 69610394, 
© 66378220, 
062704644, 
«58792290, 
034646410, 
50274967, 
« 4568 8669, 
© 49900945, 
0 3592 78655 
e30788007, 
© 25502260, 
0 2A0935 64, 
014586711, 


© 92290013, 
91031398, 
«89612366, 
» 88020507, 
«86243611, 
08427097 9, 
«62089084, 


«7816398 5, 
075381075, 
0723821615 
069143312, 
«65662770, 
261964146, 
«57981581, 
«53789863, 
06 9374505, 
«44746762, 
«39920598, 
234912599, 
229761822, 
024429591, 
018999246, 
013675835, 


041246292, 10128059, 
eO56371978, 045111106, 
DATA (TERO(ID», Tuts 76) / «18L1505S9E-08, 1 9816AG2E-08, 


© 216765 96£= 04, 
+30970439E-08, 
0 44107647E-08, 
© 6261677 2E = 08, 
© 8860 6977E-08, 
012498993E-07, 
01757 4484E-07, 
o24632041E-07, 
eo S4H13U71E-07, 
04792551 7E= 07, 
e 66529674E-97, 
©9206 0694E-07, 
0 12698235E-06, 
0174591 35E-06, 
o23928327E~ 06, 
«32689796E<06, 
04451 6637E-06, 
e60428629E=06, 
© 81766120F=- 06, 


0127927 41E-05, 
01 7171961E- 05, 
«2297 6636E- 08, 
© 30665257E- 05, 
240 762749E=05, 
©53994268E- 05, 
071327211805, 
09392324 36-05, 


02370611 8F <0 By 
»33843033E-08, 
¢48160210E-0 8) 
«64315260E-08, 
965955985 -08, 
13614673E-07, 
19127901E-07, 
e267 67A8H1E-07, 
0373 95414E=0 7, 
«320 36639E-0 7, 
7217 B923E-07, 
©99797993E=-0 7, 
01375465 8F=0 6, 
0186 96240E-06, 
«25867721 8E-06, 
o35324013E-0 6, 
o48065419E-0 6, 
+65193709E-06, 
06514322 E-06, 


e837 73933EMu 5, 
2184742S50E-05, 
e2469937KE-05, 
«329166266 -05, 
o&3727530E-05, 
«5790350 3E-0 5, 
© 76430199E-05, 
eL0056235F-04, 
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0890078126, 
0033841222, 


o25920474E 08, 
© 36974673E~08, 
e THEA THSIOE“18y 
0 10528102E -07, 
0 1882697hE-07, 
© 2081 4463F ~07, 
02912 €498F-07, 
o4062761BE-07, 
2 564869121E-07, 
e7B8292207E-07, 
e LBBLK39KE 06, 
0 t4895557E~ 06, 
«204 7548F - 06, 
0 2797 9245E- 16, 
0 38162867E~-06, 
«51 586725E-06, 
oe 70320473F-06, 
9499867 9-06, 


©1462 76 16E-15, 
0194972 328E~-05, 
e 2654596 8E-05, 
2 353492 75E <5, 
2620833535 -05, 
061681694E-05, 
010764921F 06, 


92050518, 
090760629, 
» 9304233, 
287680310, 
2 85 864995 
«63850807, 
82627102 7/ 


080676772, .801858283, 


«77620027, 
074800328, 
071753675, 
068466555, 
«64937669, 
061164127, 
057161576, 
052924362, 
© 484655%9, 
43796909, 
o 38932970, 
233890815, 
028669972, 
o 23382192, 
017901181, 
012362290, 


00786577205 . 067621596, 
0022564575, 


028 336002E-08, 
24038801 8E-08, 
oS 73OZ2144E-08, 
081266442E-1 8, 
01147264 5F -07, 
016143994E-07, 
«22645204E<07, 
03166297 7E -07, 
o44130 364KE-07, 
¢61310313€-07, 
eBRIDGBZHLE M7, 
011720776c-06, 
216128096E-06, 
o2P121 787 E06, 
eSN2KSS71E <0K, 
04122162 6E-06, 
-56N00632E-06, 
o 75835232E-06,y 
eLO236686E-05 / 


DATA (TEROCTID, 19795156) / o1L0Z8445E-05, .11879068F-05, 


015958281E-05, 
02136978 2E-0 5, 
» 2852491 HE -05, 
oe STIS AILSE-05, 
oSHI3BBATE-05, 
066551620E-05, 
0 87704910F-05, 
eLI5S2L246F 04, 


0021283416, 0.4/7 


118 


120 


125 


130 


135 


140 


145 


150 


160 


165 


POONA UFUNY DO 


0 12326243E-04, 
016130192E-06, 
e210372Z10E= 04, 
027349351E- 06, 
e S544 18SLE- Ob, 
«457862082E- 04, 
0586980187E= 04, 
oe 7566. 3267E= 0b, 
© 9680 &43KE- 06, 
012345698E-03, 
0 15694459E-03, 


0131 6912 7E=8 &, 
01724276 BE-Dh, 
022470265" -04) 
o29189025F-9 4, 
03779568 5E-04, 
o& 876353 2E=0 4, 
0S2764699E-04, 
o 8049481 7E“0 by 
¢10290362E-03, 
0131130%7E-03, 
¢ 1665661 9F-03, 


e 141073065 -64, 
0 18428397E= 06, 
0239962 JSE-04, 
o SLL4EZ17E= Oh, 
240297636E ~06, 
251971360E- 0b, 
0668126745 - 04, 
2 85617765E~ 0h, 
0109365216 = 03, 
0139253 05E-03,5 
017674231E-03,5 


015086387E-06, 
01969261 SE-04, 
02562050 0F 04, 
033227997E=04, 
eb2986707E-fh, 
0$5356629E-04, 
072107499E-04, 
69104854 2E-04, 
011620929E-03, 
© 1478491 9E-03, 
018750262E-03 / 


DATA (TERD(L)» T#157,234) /¢ .19887824E-03, «210901 84E-03, 


PORNO EFUNKrPOBanGwRFeUne 


0 2236 0761E- 03, 
0 28199597E= 03, 
6 3546 9649E-03, 
0 44420794E~ 03, 
2 5848 4722E- 03, 
069082932E-03, 
«657386992E-03, 
010607098E-02, 
1308 0500E-02, 
©16079137E-02, 
14979 2048E~ 02, 
© 240641365-02, 
©29298098E=02, 
2 35556487E-02, 
o43013869E~02, 
«528690676 = 02, 
06234744 0F-02, 
747031 768-02, 
6 89221551F-02, 


e237 03144E 83, 
©2986 85 96E=-0 2, 
©37517508E-03, 
e46974632E-0 3, 
o58627725E-8 3, 
o72937850E-03, 
~ 904 509%9E-0 3, 
011186107S5E-02, 
013777304E-02, 
01692212 366-02, 
020716409E-8 2, 
025285285" 02, 
o30760230F-02, 
0 3730109ZE-0 2, 
0450 86292E-0 2, 
»5432706F -02, 
e6S52K9776F-0 2, 
o 76118166E-02, 
093225623E-02, 


2 25121089E-03, 
+ 316300S3E-03, 
039698 27HE 08, 
°49665360E-03, 
© 61936378E-03, 
2769926 76F-03, 
0 95402778E-83,5 
oe 11783764E-02, 
oL GS5OR32Z5E -02, 
01 7B80S77IE “02,5 
oe 217828462E~-02, 
o 2656 308% —- 02, 
o S2288671E -02, 
e 391236 73E 02, 
» 47253306E-02 5 
e 56890172E-02, 
06827 3562E=-82, 
oe 816729306 -02, 
o 97369910502, 


«2661853 3° -03, 
o $34 8868 6F-03, 
41997 4006-03, 
25249971 3E-03, 
«6541 6671E-03, 
06125624 7E-03, 
© 10060558E-02, 
0 12616455E-02, 
015278078E-02, 
018771 G00E~02, 
oe 2289738 3F-02, 
027899886E-02, 
2 33886700F-02, 
«41026681E-02, 
o&951237 46-02, 
¢ 5956228 7E=-02, 
©71423190E-02, 
08537237 8E-02, 
010171986E-01 / 


DATA (TERDCI) » 182359312) / .L1062210GE-02, 1108993 0E=01, 


GAN HPO SNFUFUNE DANFREUN YP 


01157 6041E-01, 
01371 56S0E- 01, 
0161986606E- 01, 
¢19070802E-01, 
0 2237932%E- 01, 
o 2617647S5E~-01, 
o 30524740E-01, 
o35478877E~Diy 
0 41105318E-01, 
047471879E- 01, 
0$4649361E-01, 
062711080E-01, 
0 7173204CE-01, 
©61780571E-01, 
09295 7046E-01, 
010531307E +00, 
01189 3029E+ 08, 
013367915E+00, 
015022550E+ 00, 


012081 0G3E-01, 
oAH302Z545F 01, 
1168 764465 -61, 
»19854636t-01, 
023280999E-01, 
o 27211441 F-81, 
03170386316 -81, 
036819865E-01, 
04262485 SE -81, 
0491 87401E-01, 
565 78979E-01, 
0646733 98E-01, 
07616620 3E-01, 
0846 73470E-01, 
095931799F-815 
© 10859632: +00, 
122540006400, 
01376322 7E+00, 
015433 761E +00, 


2126055546F-01, 
0149115716 -01, 
2 175836 88F-01, 
o 2066698SE-01, 
024216156F-01, 
0 28279S10E=-01,5 
0 32921 6B1E- 01, 
e 38203897F-61, 
046191723F -81, 
= 50954726E-01, 
» 38565 016€~01, 
©67096876E- 01, 
o 766260 82E-Gi» 
o 8722 90S9E-01, 
0969819S1E-01, 
0 11195954E +00, 
© 126236 02E +88, 
0161373 7AE 408, 
0196941715408, 


013150207E-01,5 
o15543422E-01, 
016323482E-01, 
©21508270E-01, 
0251 79665E-01, 
0293836246 -01, 
0 3417989ZE-01, 
e 39632 028E-01, 
eb SBN7OZ7E“O02,5 
052774996E-01, 
26 0608646E-31, 
06 93862697E-01, 
e79173173E-01, 
2 90056524E-01, 
010219666 +00, 
011540 363E+00, 
0130023606400, 
01 &600651E+00, 
016 363862E 000 + 


DATA (TERDCI)» te 313,390) / .16802916E+08, .17272408E+ 00, 


027749413E+ 00, 016236999E 400, ~187SH2SZECOD, 01924117 SE 400, 


019757680E +60, 
0219233536 400, 


«202844062 +60, 
e224 89567E+0 0, 
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o 2062 O799E +00, 
o 2306 S833E +00, 


©21367101E +00, 
o23652112E 906, 
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2 Ee ee EOE 


170 


175 


180 


185 


490 


195 


210 


215 


40 


PUSAN EFP COIOVE UNS SFPBWOEOVNAUSUN- CAVNFuTs 


ERF 
ARG 


24248683E+ 00, 
+26734435E 000, 
293811 2%E+ 00, 
03218 G671E +005 
+ 3514 7662 +085 
«3625799064005 
41510750808, 
«468961678 608, 
oh BbO2506E +00, 
«S201 6058E+00, 
58720760600, 
+59498579£+68, 
6 332957KE 6 06, 
067191 881E+ 00, 
«7 1061978E +005 
749148726600, 


08060 433KE +00, 
© 84297518E +00, 
© 8787 62S8E +08, 
091318411E¢00, 
© 94590063E +00, 
«97665698 +00, 
¢10052957E+01, 
©1031 2609E¢ 01, 
0 10546218E+ 01, 
«107506626401, 
0 10924056E+ 02, 
¢ 11064763E+01, 
04217151266 ¢ 01, 
© 12263263E +01, 
o11279279€+¢ 01, 
= 1. 


MX 2 &.%% 


GG TO 10 
ENTRY ERFD 


ERF 


ARGHX = &.5 


vos 


Ags(%) 


SIGN = 1. 
TFCRLT.0.) SIGN = =f, 


If 
= 


(ARGHX aL Fe 
100.°V + .5 


M2 © .O3°FLOAT CL) 
OX = VoxXt 


Jj 2 


&$1-L 


0248546365260, 
02736108 4E+80, 
¢ 30067728E+0 0, 
o32912497E+8 G, 
0339113156900, 
039056184E+00, 
o4234506SE+6 6, 
045761 $25E 000, 
ob929EH67F 400, 
0352936346 BE+0 0, 
2 36659089- +00, 
06045216 7E*0 0, 
06429310 7E*00, 
068159626" 60 0, 
© 7202761 9E€+00, 
07387235 8E+00, 


o81536635E+00, 
285204344 +8 C, 
0887S2573E+0 0, 
292153201: +00, 
09337642 7E+00, 
© 964011436400, 
0102195316901, 
01037 3633F+01, 
0106 001%1E+01, 
2 107 969896 +01, 
01096235 7E+01, 
0 110946985 +0 2, 
011192 76 26+0 1, 
o12255617E+04, 
0112 82663E+01, 


Vd RETURN 


IF CERF) 20, 30, 20 
ERF s(TERACJS-36) © TERD(J) SOX* (2.°0X"X1)) SSEGN 
RETURN 
ERF S(TERDC I) © (Le = OX*(DX. & 2o*K15(2e-OK9KL) ) DI*SIGN 


RET 
ENO 


URN 
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o2547130ZE +005 
©28037767E +00, 
0 30764223E +80, 
0336479606 +00, 
o 3668 6% 30E +00, 
0 39667240E+00, 
o43187S557E +60, 
2466350516 +00, 
-5U196657E 408, 
o$3854079E+00, 
© 57601720F +08, 
0614067566400, 
©6525 6247E+08, 
-691276£86E°00, 
o72992166E +00, 
0 76826714E +00, 


0 8246 3224E+ 005 
2961037 03F +00, 
o 896176626 +00, 
¢92977025E +00, 
0 96154130F 400, 
¢ 991221006 +00, 
of0165522E+01, 
¢ 10432932E 401, 
0 10652209E +01, 
0 10861388E +01, 
018998593E+01, 
¢ 11122469401, 
9 11211606E 401, 
¢ 112657526 +01, 
01128 3792E+01 / 


«26097837E+00, 
o Z26706457E+00, 
o 314705570+00, 
e 3439297 GE +00, 
0 37466696E 400, 
o 40684712 3E+00, 
eo 4603797 9F +00, 
o47515331E +00, 
0511034716400, 
«5478697 2E +00, 
e 5$8548322E +60, 
262368006E 400, 
» 662246286400, 
»70095067E+00, 
0 7395467 6E 200, 
o77777518E¢00 / 


DATA CTERO(I)» 183919451) / o7B724343E+00, « 79666759E400, 


e 83383 665E+00, 
o 869951555 200, 
©90073687E 400, 
o 937A94USKE 000, 
© 96916756E+00, 
0 998283725 280, 
010249696E+61, 
0 1049047 26201, 
010702393E001, 
© 10883712€+01, 
0110327H1E 004, 
0411480816401, 
0112206368401, 
0112736412901, 


SUBROUTINE PSXY (XoVoIXsIY¥ ,ANS) 
COMMON /T101/U,V5 PHI» CEE »SIGMAXs SIGMAY REI LL» SNPHI,CSPHE 
QzyY 
T1izCOS(xX) 
5 T3aSIN(X) 
AZUeR*TL*CSPHI-CEE*R* T3*SNPHI 
Bev+R T1i*SNPHT +C EE*R® T3*CSPHI 
AZA5S CA) 
B= ASS (8) 
10 CALL FKILL ( RKILL»SIGMAXsSIGMAY, Ay, ANS ) q 
ANS#®* ANS 
RETURN 
END 


SUBROUTINE PBXYS(XsV¥pIXeIV,ANS) 
COMMON /T101/U,V, PHT, CEE, SIGMAX, SIGMAY »RKILLy SNPHI pCSPHI 
91 FORMAT € 1H095H IR=yI7, 29H CIRCV HAS BAO INPUT VALUES U 
RsY 
5 Tie COS(X) 
T32SINCX) 
ASU+R*T 1°CSPH I-CEEFR® T3°SNPHI 
B2V¥ +R*T1* SNPHIFCEE*R®* 13°C SPHI 
O4=SQRT (¢ AT ACB*B ) 


10 CALL CIRCV ( RKILL»O1,t,ANS,IR ) 
TF € IRNE.O > PRINT 92,fR 
ANS sANS*R 
RETURN 
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